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ABSTRACT 


The  Matrix  Pencil  Approach  [1]  was  shown  to  be  an  effective  and 
efficient  method  for  estimating  the  angles  of  arrival  of  multiple  narrow- 
band  sources.  It  is  classified  as  a  non  search  procedure.  Therefore,  it  is 
computationally  less  complicated  and  eliminates  problems  encountered  in 
search  procedures  with  regard  to  memory  storage  and  system  calibration. 
Having  collected  the  data  from  the  outputs  of  a  linear  uniformly  spaced  ar¬ 
ray  consisting  of  m  sensors,  the  objective  is  to  estimate  the  locations  of 
the  d  sources  (d<m).  The  information  about  the  parameters  of  interest  are 
contained  in  the  rank  reducing  values  of  a  matrix  pencil  generated  from  the 
set  of  data. 

Several  extensions  of  the  Matrix  Pencil  Approach  appear  in  this 
work.  In  the  earlier  work  [1],  a  data  window  of  length  L=  m-d  was 

used  to  form  (d+1)  vectors.  Because  this  choice  results  in  a  minimum  number 
of  vectors  to  span  the  array,  it  fails  to  take  into  consideration  the  pos¬ 
sible  separation  of  the  signal  and  noise  subspaces.  Thus,  performance  is 
drastically  degraded  at  low  values  of  signal  to  noise  ratio  (SNR).  In  this 
work  it  is  shown  that  improved  performance  can  be  achieved  using  a 

data  window  of  length  L=d.  Because  this  results  in  (m-d+1)  vectors,  which 
is  the  maximum  number  of  vectors  to  span  the  array,  identification  of  the 
noise  subspace  is  possible.  Previous  developments  of  high  resolution  algo¬ 
rithms  neglected  the  effects  of  mutual  coupling  which  occurs  between  the 
elements  of  an  array.  We  show  that  failure  to  account  for  mutual  coupling 
results  in  poor  performance.  Concentrating  our  efforts  on  the  non  search 
procedures  of  ESPRIT  and  the  Moving  Window,  we  have  successfully  improved 
the  performance  of  these  methods  by  compensating  for  the  mutual  coupling. 
The  problem  ,f  wideband  signals  is  much  more  difficult  and  has  been  studied 


i 


by  only  a  few  investigators.  Three  methods  dealing  with  the  wideband  case 
are  proposed  in  this  work.  ,  In  the  first  method  the  wideband  signals 

are  modeled  as  sums  of  decaying  exponentials.  This  model  is  suitable  for 
non  stationary  signals.  The  natural  frequencies  of  the  sources  are  assumed 
to  be  unknown  at  the  receiver.  Therefore,  the  estimation  procedure  consists 
of  estimating  both  natural  frequencies  and  angles  of  arrival  of  the  sources 
by  means  of  two  matrix  pencils.  However,  an  ambiguity  problem  arises  as  to 
which  natural  frequencies  are  to  be  associated  with  which  angles  of  ar¬ 
rival.  We  show  that  a  third  matrix  pencil  removes  this  ambiguity.  A  second 
method  is  proposed  where  the  wideband  sources  are  assumed  to  be  linear  sys¬ 
tems  driven  by  white  noise.  This  model  is  appropriate  for  stationary  sig¬ 
nals.  The  same  array  configuration  as  in  the  first  case  is  used.  The  analy¬ 
sis  is  carried  out  on  the  unit  circle  using  the  Discrete  Fourier  Transform. 
The  third  approach  makes  use  of  the  coherent  signal  subspace  method  (CSS) 
proposed  by  Wang  and  Kaveh  [56]  in  conjunction  with  the  moving  window  oper¬ 
ator.  Again  the  method  performs  relatively  well  when  compared  to  ESPRIT. 
Finally,  we  have  studied  the  effects  of  perturbation  due  to  noise  and  due 
to  sensor  spacing.  We  have  derived  upper  bounds  for  the  Chordal  Metric 
which  is  a  measure  between  the  true  eigenvalue  and  the  perturbed  one.  The 
chordal  metric  is  shown  to  be  a  functional  of  the  true  and  the  perturbed 
angles  of  arrival. 

Computer  simulations  are  carried  out  for  each  of  the  analyses  as¬ 
sociated  with  respect  to  data  window  length,  mutual  coupling  compensation, 
the  three  wideband  methods,  and  the  upper  bounds  in  the  chordal  metric. 
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CHAPTER  1 


INTRODUCTION 


1.1  RESEARCH  OBJECTIVES 

In  this  work,  the  problem  of  passive  high  resolution  direc¬ 

tion  finding  of  multiple  sources  using  a  sensor  array  is  addressed.  This 
problem  arises  is  such  systems  as  radar,  sonar,  seismology,  geophysics, 
etc.  A  direction  finding  system  is  referred  to  as  passive  when  the  signals 
received  at  the  array  are  generated  externally  to  the  array.  These  signals 
can  be  either  narrowband  or  broadband.  In  both  cases,  given  measurements 
collected  at  the  array  output,  the  objective  is  to  determine  the  number  of 
targets  (Detection)  and  estimate  their  parameters  such  as  angles  of  ar¬ 
rival,  natural  frequencies,  etc.,  (Estimation).  A  signal  is  classified  as 
narrowband  when  the  bandwidth  of  the  impinging  signals  from  the  sources  is 
much  less  than  the  reciprocal  of  the  propagation  time  of  the  wavefronts 
across  the  array.  When  this  condition  does  not  hold,  the  signals  are  said 
to  be  wideband. 

The  problem  of  narrowband  sources  has  been  studied  «  tensively. 
Solutions  range  from  the  classical  ones  such  as  the  periodogram,  the  cor- 
relogram,  etc.  to  subspace  approaches  such  as  MUSIC,  ESPRIT,  Matrix  Pencil, 
etc.  This  work  deals  with  the  Matrix  Pencil  Approach  [1].  For  the  sake  of 
clarity,  assume  that  d  narrowband  sources  are  present,  m  measurements  are 
collected  at  the  output  of  a  linear  uniformly  spaced  array  of  m  elements 
and  m>d  (Fig.  1-1).  The  Matrix  Pencil  approach  is  based  on  an  invariance 
introduced  by  the  geometry  of  the  array  (linear  uniformly  spaced).  It  is 
shown  that  the  parameters  of  interest,  i.e,  the  angles  of  arrival,  are  re¬ 
lated  to  the  rank  reducing  values  of  a  matrix  pencil  generated  from  the 
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data.  The  method  is  shown  to  hold  even  in  the  presence  of  fully  correlated 
sources.  However,  the  method,  as  applied  in  [1],  did  not  fully  take  into 
account  the  separation  of  the  noise  and  signal  subspaces.  The  eigenvalues 
generated  from  the  data,  as  applied  in  [1],  were  obtained  using  a 
dimensionality  too  small  to  allow  for  identification  of  the  noise  subspace. 
That  is  the  reason  why  the  method  performed  relatively  poorly  at  low  signal 
to  noise  ratios  (SNR).  We  have  devised  a  scheme  which  displays  improved 
performance  at  low  SNR.  High  resolution  algorithms  devised  previously  ig¬ 
nored  mutual  coupling  which  may  exist  between  array  elements.  In  effect, 
each  element  of  the  array  was  modeled  as  though  it  existed  by  itself.  Ue 
have  studied  the  effects  of  mutual  coupling  and  we  devised  schemes  to  com¬ 
pensate  for  these  effects.  The  case  of  wideband  sources  is  much  more  com¬ 
plex  and  has  been  treated  by  only  a  few  authors.  In  this  dissertation,  we 
extended  the  notion  of  the  matrix  pencil  to  the  wideband  case.  Finally, 
previous  algorithms  assumed  a  perfect  sensor  spacing.  We  studied  the  per¬ 
formance  of  the  matrix  pencil  for  the  case  of  small  perturbations  in  the 
sensor  spacing. 

1.2  SIGNAL  MODEL 

Before  reviewing  the  work  that  has  been  done  previously,  it  is 
useful  to  develop  expressions  for  the  received  signals  and  see  how  they 
simplify  for  the  case  of  narrowband  signals.  For  this,  assume  that  we  have 
a  linear  uniformly  spaced  array  composed  of  m  identical  sensors.  Let  A  be 
the  sensor  spacing  and  d  the  number  of  sources.  These  sources  are  assumed 
to  be  in  the  far  field  so  that  planar  waves  arrive  at  the  array.  It  follows 
that  the  output  of  the  i-th  array  element  can  be  expressed  as 
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(1.2-1) 


d 

x^(t)  ■  T  ^( ®k)sk( t— x^k)  +  n^( t)  j  i  =  l»  2,  -  •  .,  m, 

k-1 


where 

a(©k)  is  the  gain  pattern  of  the  sensor  at  angle  0^, 
nj(t)  is  the  additive  noise, 

Xjk  is  the  time  delay  that  source  k  takes  to  travel  from  the 
reference  sensor  to  the  i-th  sensor.  With  respect  to  the  first  sensor  x^ 
is  given  by 

xik  *  (i-l)(A/c)sin(0k),  (1.2-2) 

where  c  is  the  propagation  speed  of  the  plane  waves. 

Let  sk(t)  be  a  modulated  signal  of  the  form 
sk(t)  -  gk( Ocos^t+a^t)),  (1.2-3) 

where  all  the  sources  are  assumed  to  be  emitting  at  the  same  carrier  fre¬ 
quency  uq.  Therefore,  sk(t-tjk)  is  given  by 

sk<  t-Tik)“Sk(  t-Tik)cos(«ot-«o‘Cik+<*k(  t~Tik)  >  •  (1-2-4) 

Note  that 

«^Tik*(i-l)^k,  (1.2-5) 

where 

♦k  -  «q( A/c)sin(0k) .  (1.2-6) 

For  narrowband  signals,  the  modulation  varies  slowly  relative  to  the  car¬ 
rier.  In  particular,  assume  that  gk(t)  and  ^(t)  are  essentially  unchanged 
over  the  duration  of  the  observation  interval.  Then 

«k<t“Tik)  s  Sk^)  (1.2-7) 

and 

®k( t-Tik)  a  «k(t).  (1.2-8) 

The  expression  for  Sjt(t-Xik)  simplifies  to 

sk(  t_Tik)*«k^  t)cos(«ot-(i-l)4»k+ak(  t)).  (1.2-9) 
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Clearly,  for  narrowband  signals,  we  see  that  the  time  delay  results  in  a 
phase  shift.  When  the  signals  sk(t)  are  broadband,  gk(t-Tik)  and  ok(t-xik) 
cau  no  longer  be  approximated  by  gjc( t )  and  otk(t),  respectively. 

For  narrowband  signals,  note  that  s^t-x^)  can  be  written  as 

sk<t-‘cik>  “  Ret  Sk^)  eJ^*  }  (1.2-10) 

where  Re{.}  denotes  the  real  part  operator.  Let  sk(t)  be  the  complex  en¬ 
velope  of  s^t).  Then 

sk(t)  -  gkCtJeJ0^^).  (1.2-11) 

Also  define  xj(t)  and  n^(t)  as  the  complex  envelopes  of  x^(t)  and 
nj(t),  respectively.  It  is  easy  to  see  that 

d 

Xi(t)*  £  ar  \)sk_(t)e-J<1-1>^k  +  n4(t)  ;  i-1,  2,  .  .  . ,  m.  (1.2-12) 
k=l 

In  the  remainder  of  this  dissertation,  we  will  drop  the  and  will  assume 
that  we  are  dealing  with  the  complex  envelopes  of  the  corresponding  sig¬ 
nals.  The  expression  for  the  received  signal  at  the  i-th  sensor  will  be  ex¬ 
pressed  as 

d 

Xi(t)  -  £  a(ek)sk(t)e-j<1-1)*k  +  n^t)  ;  i-1,  2,  .  .  .,  m.  (1.2-13) 

k-1 

1.3  LITERATURE  SURVEY 

The  problem  of  estimating  the  angular  locations  of  sources  is  of 
great  importance  and,  over  the  years,  has  occupied  many  researchers.  This 
problem  is  the  spatial  frequency  analog  of  the  temporal  problem  dealing 
with  harmonic  retrieval  in  additive  noise.  Consequently,  research  in  direc¬ 
tion  finding  has  benefited  a  great  deal  from  the  advances  made  in  spectral 
analysis.  The  periodogram  [2,3]  was  seen  as  one  of  the  most  promising  meth- 
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ods  in  determining  the  locations  of  sources  where  one  has  to  plot  the  func¬ 
tion 

d 

Pxx(*)  -  (d/m)  |  E  xi(t)e-J1+  |2,  (1.3-1) 

k>l 

where  xj(t)  is  the  received  signal  at  the  i-th  sensor,  A  is  the  sensor 
spacing  and  m  is  the  total  number  of  sensors.  The  periodogram  has  the  ad¬ 
vantage  of  being  non  parametric  in  the  sense  that  it  does  not  rely  on 
knowledge  of  a  model  of  the  input  processes.  Also,  it  is  robust  in  that  it 
is  relatively  insensitive  to  signal  parameters.  It  is  also  simple  to  imple¬ 
ment.  However,  a  disadvantage  is  that  the  physical  size  of  the  array  has  to 
be  increased  in  order  to  achieve  a  better  resolution.  Typically,  two 
sources  that  are  separated  by  less  than  one  standard  beamwidth  cannot  be 
resolved.  The  standard  beamwidth  is  defined  as  ♦q  =  2n/m. 

Another  way  of  estimating  the  power  spectral  density  (p.s.d)  is  by 
using  the  autocorrelation  sequence  [4,5]  defined  as 

(m-1) 

rxx(l)  =  1/ (m-1)  Z  xn+1(t)xn*(t)  ;  1=0,  1,  .  .  .,  (m-1).  (1.3-2) 

n=l 

When  considering  only  a  finite  sequence,  the  power  spectral  density  is 
estimated  by  the  correlogram  which  is  given  by 

L 

Pxx(*)  =  A  Z  rxxaU~il*  ,  (1.3-3) 

1— L 

where  L  *  (m/10)  has  been  found  to  give  good  estimates  of  the  terms  in  the 
autocorrelation  sequence.  Note  that  this  may  require  an  unacceptably  large 
value  of  m.  This  value  of  L  arises  in  an  attempt  to  get  good  estimates  of 
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rxx(l)  for  all  lags  in  the  sum. 

Parametric  spectral  estimation  techniques  [6,7]  have  been  intro¬ 
duced  to  overcome  the  limitations  encountered  with  the  periodogram  or  the 
correlogram.  Here,  the  spatial  samples  x^(t)  are  taken  as  the  samples  from 
an  autoregressive  (AR)  process  of  order  p  which,  by  definition,  satisfies 
the  linear  difference  equation 

P 

xn(t)  =  -  Z  akxn_k(t)  +  un  ,  (1.3-4) 

k-1 


where  the  coefficients  ak's  ;  k=l,.  .  .,p,  are  constant  parameters  and  un 
is  a  sample  from  a  zero  mean  white  Gaussian  process  with  autocorrelation 
sequence 

r  E[ |un|2J-Puu  ;  1-0 

ruua)«  (1.3-5) 

{  0  j  1*0 


In  this  approach,  the  angles  of  arrival  are  obtained  as  the  maxima  of  the 
power  spectral  density 


PAR<*>  -  (APuu> 


1+ 


P 

Z 

k-1 


ake-3k* 


(1.3-6) 


where  puu  and  ak  ;  k-1, 2,.  .  .,p  are  obtained  by  solving  the  Yule-Walker 
normal  equations 


rxx(0)  rxx  (1)  .  .  .  rxx  (p) 

■  1  ■ 

Puu 

rxx(*)  rxx(®)  •  •  •  rxx  (P~l) 

•  «  «  •  •  • 

al 

= 

0 

•  •  •  •  •  • 

■  rxx(p)  rXx(P“l)  •  *  *  rxx(0) 

■ap  ■ 

.  0  . 
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and  rxx(l)  is  the  autocorrelation  estimator  given  by 


(m-1) 

rxx(l)  «  1/ (m-1)  E  xn+1(t)xn*(t)  ;  1=0,  1 . p.  (1.3-7) 

n=l 


Recall  that  the  p.s.d  in  the  correlogram  method  is  given  by 

L 

Pxx(*)  ,  A  E  r^CDe-J1*  ,  (1.3-8) 

1—L 


which  means  that  the  autocorrelation  sequence  is  assumed  to  be  zero  for 
|l|>L.  This  windowing  is  the  reason  for  the  poor  resolution  capability  of 
the  classical  estimator.  The  AR  p.s.d  estimator  uses  the  autocorrelation  of 
the  correlogram  methods  and  extrapolates  estimates  of  the  autocorrelation 
sequence  through 


L 

rxx(l)  *  E  a^r  xx(l-k)  ;  1>0 
k»-L 

rxx(l)  *  r  XX^l)  » 
which  results  in  a  p.s.d  given  by 


(1.3-9) 


Pxx(*)  ,  A  E  r^l^1*  .  (1.3-10) 

l»-» 

This  means  that  it  extrapolates  estimates  of  the  entire  autocorrelation  se¬ 
quence  which  explains  the  high  resolution  property  of  the  AR  p.s.d 
estimator.  Burg  [8]  showed  that  the  extrapolated  autocorrelation  function 
has  maximum  entropy.  This  results  in  the  most  random  time  series  which  has 
rXx(0)>  rxx(l)*  •  •  >i  rxx(p)  as  *ts  first  (p+1)  lag  values. 

Linear  prediction  techniques  [9-11]  can  also  be  applied  to  the 
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direction  finding  problem.  The  idea  here  is  to  estimate  the  output  of  the 
n-th  sensor  as  a  linear  combination  of  the  other  sensor  outputs  ;  i.e, 

L 

xn(t)  •  -  I  afk  xn_k(t)  5  L<n<m  (1.3-11) 

k=l 

where  L  is  the  order  of  the  prediction  filter.  The  coefficients  a^k  are 
chosen  such  that  the  error  p=E[ |xn(t)-xn(t) |2]  is  minimized.  The  mini¬ 
mization  results  in  the  minimum  error  variance  given  by 
Pmin  3  rxx<°)  +  ELHa- 

The  equations  involved  in  the  minimization  can  be  expressed  as 


rxx(0)  rxx  ( 

fxx*1)  ^xx<°) 

•  • 

.  .  .  rxx*(L)  ‘ 

.  .  .  rxx*(L-l) 

•  •  •  • 

'  1  ' 

4 

• 

H 

PUU 

0 

• 

•  • 

•  • 

rXx<L>  rxx^-D 

•  •  •  • 

•  •  •  • 

.  .  .  rxx(0) 

• 

af 

L  »p  J 

• 

• 

.  o  . 

These  equations  are  identical  to  the  Yule-Walker  equations  when  L=p.  There¬ 
fore,  for  L-p,  the  Yule-Walker  equations  arise  in  both  the  AR  and  linear 
forward  prediction  approaches.  Backward  linear  prediction  can  also  be  in¬ 
troduced  where 

L 

Xn-L<l>  ■  *  1  abk  Xn-L-k^)  •  (1.3-12) 

k-1 

Again,  the  coefficients  a^k  are  chosen  so  as  to  minimize  the  error  defined 
as  p-E[ |xn_L(t)-xn_L<t) |2] .  This  leads  us  to  the  same  set  of  equations  as 
before  and  it  can  be  shown  that  P^min  =  P^min  an<*  (a^k)=(abk)*  ^or  ^=1,  2, 

.  .  .,  L.  The  forward-backward  linear  prediction  (FBLP)  is  based  on  the  use 
of  least  squares  for  estimating  the  AR  parameters.  This  technique  is  also 
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known  as  the  least  squares  method  [12-15].  Here,  ve  have  to  solve  an  over¬ 
determined  set  of  equations  of  2(m-p)  equations  in  p  unknowns,  where  p  is 
the  order  of  the  filter. 

If  the  available  data  sequence  is  very  long,  sequential  estimation 
techniques  are  available  for  updating  the  AR  parameter  estimates.  These 
techniques  are  useful  for  tracking  sources  with  slowly  varying  angles  of 
arrival.  They  are  also  known  as  adaptive  algorithms.  The  least  mean  squares 
(LMS)  adaptive  algorithm  [16]  is  based  upon  the  gradient  steepest  descent 
adaptive  procedure  where  the  (j+l)-th  element  is  given  by 

P 

Xj+j(t)  *  -  I  ak  xj+^_k(t)  .  (1.3-13) 

k=.l 

A  minimization  procedure  is  then  developed  which  results  in  a  (j+l)-th  er¬ 
ror  given  by 

ej+1(p,j)-xj+1(t)+(xj(p-l))Tap(j).  (1.3-14) 

The  LMS  algorithm  attempts  to  find  the  minimum  of  the  mean-squared  error 

quadratic  surface.  This  search  proceeds  in  a  random  fashion  but,  on  the 
average,  the  LMS  algorithm  converges  to  the  optimum  coefficient  vector.  The 
condition  for  convergence  is  that  the  step  size  satisfies 

0<*«l/prxx(0)),  (1.3-15) 

where 

rxx(0)  «(l/j)  Z  |xk(t)|2  .  (1.3-16) 

k-1 

Another  way  of  dealing  with  long  data  sequences  is  to  use  a  recur¬ 
sive  least  squares  method  (RLS)  [17].  This  method  searches  for  the  forward 
prediction  error  filter  vector  ap  which  minimizes  the  sum  of  the  squared 
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errors  subject  to  the  constraint  that  the  first  component  of  is  unity. 
This  method  leads  to  the  same  answer  as  in  the  LMS  case.  The  only  dif¬ 
ference  is  that  the  adaptive  gain  is  constant  for  the  LMS  whereas  it  is  a 
spatially  variant  matrix  for  the  RLS.  The  LMS  appears  to  be  more  attractive 
than  the  RLS  because  of  its  robustness  in  its  behavior,  its  insensitivity 
to  perturbations  and  its  computational  flexibility.  It  requires  a  number  of 
computations  proportional  to  p  whereas  a  number  of  computations  propor¬ 
tional  to  p^  is  required  for  the  RLS. 

We  have  seen  earlier  how  an  AR  process  can  be  used  to  generate  ob¬ 
served  data  samples.  The  spatial  samples  x^(t)  can  also  be  modeled  as 
though  they  are  generated  by  a  moving  average  (MA)  process  [18]  where  by 
definition,  we  have 


q  q 

xn(t)  *  Z  bk  un_k  +  un  =  I  bk  un_k  ,  (1.3-17) 

k=l  k=0 


where  bg=l,  q  is  the  order  of  the  MA  process  and  un  are  samples  from  a  zero 
mean  white  noise  process  with  autocorrelation  sequence  given  by 


fuud)- 


E I  I  un  I  ^ 1 =  Puu 
0  ;  1*0 


1=0 


The  power  spectrum  is  then  given  by 


?(i)  «  (AA^) 


1+  E  bke-jk* 
k=l 


(1.3-18) 


where  ♦*(<oD/c)sin(0) .  The  angles  of  arrival  are  determined  as  the  peaks  of 
this  spectrum.  To  avoid  solution  of  a  large  number  of  nonlinear  equations, 
the  parameters  bk  ;  k=l,  2,  .  .  . ,  q,  can  be  estimated  by  approximating  the 


11 


MA  process  with  an  equivalent  AR  process  of  order  p  where  p»q.  The  parame¬ 
ters  of  these  two  processes  are  related  by 


q  f  1  5  1=0 

ai  ♦  Z  bn  $(1)  =  <  (1.3-19) 

n-1  I  0  ;  1/0. 

It  can  be  shown  that  the  MA  model  is  not  a  high  resolution  spectral 
estimator  because  it  does  not  model  narrowband  spectra  very  well. 

AR  and  MA  processes  can  be  combined  to  form  an  ARMA  process  [20]. 
Here  the  spatial  samples  xn(t)  are  modeled  as  samples  from  the  process 


P  q 

xn(t>  *  ~  ak  xn-k  +  ^  ^k  un-k  *  (1.3-20) 

k-l  k=0 


where  bg-1,  aj_,  a2,  .  .  . ,  ap  and  bj,  b2»*  .  . ,  bq  are  constant  parameters. 
un  is  a  sample  from  a  zero  mean  white  Gaussian  process  with 


f  E[KI2]»Puu  ?  1-0 

rUU^)*  1 

l  0  ;  1/0 


The  power  spectrum  is  given  by 


q 

i+  z 

k-l 


P 

1+  Z 
k-l 


ake 


-jk$ 


(1.3-21) 


The  angles  of  arrival  of  the  sources  are  determined  as  the  peaks  of  this 
spectrum.  Again,  to  avoid  having  to  solve  a  set  of  nonlinear  equations,  the 
Coefficients  aj,  a2»  •  •  . ,  ap  and  bj_,  b2»  •  •  . ,  bq  can  be  estimated  by 
approximating  the  ARMA  process  with  an  equivalent  AR  process  of  order 
r»(p+q).  A  least  squares  procedure  can  then  be  used  to  solve  for  the 
moving  average  coefficients  1,  b^q),  b2(q),  .  .  .,  bq(q).  These  ele- 
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ments  are  then  used  to  solve  for  the  coefficient  1,  a^(p),  a2<p),  .  . 

. ,  ap(p)  knowing  that 

q 

al  *  C1  +  ^  ^n  cl-n  (1.3-22) 

n*l 

where  c^  are  the  parameters  of  the  equivalent  AR  process. 

The  minimum  variance  spectral  estimator  (MVSE),  also  known  as 
Capon's  maximum  likelihood  estimator  [25],  does  not  make  use  of  the  stan¬ 
dard  maximum  likelihood  estimate  (MLE).  Instead,  a  constrained  optimization 
problem  is  solved.  The  MVSE  generates  a  spectral  estimate  that  describes 
relative  component  strengths  over  spatial  frequency,  but  it  is  not  a  true 
p.s.d  estimator.  Let  be  the  response  of  a  transversal  filter  with  input 
Xi(t).  yj  is  given  by 

P 

yi  -  E  an  xj.^t)  , 
k»0 

where  a^  ;  k»l,  2,  .  .  . ,  p,  are  the  filter  coefficients.  These  coeffi¬ 
cients  are  selected  so  as  to  minimize  the  variance  p=E[|yj|^],  subject  to 
the  constraint  that  a  desired  spatial  sinusoidal  input  does  not  experience 
distortion.  It  should  be  pointed  out  that  this  estimator  is  a  non 
parametric  estimator. 

The  beam  forming  (BF)  [27]  algorithm  is  suited  to  single  sources 
where  in  effect,  the  observations  are  modeled  as  an  MA  (Moving  Average) 
process. 

Recently,  subspace  approaches  have  been  introduced.  These  methods 
are  based  on  an  eigenvalue-eigenvector  decomposition  of  the  correlation 
matrix.  This  makes  use  of  the  fact  that  there  is  a  relationship  between  the 
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eigenvectors  of  the  spatial  correlation  matrix  and  the  source  angles  of  ar¬ 
rival.  These  methods  are  shown  to  yield  asymptotically  unbiased  estimates 
of  such  signal  parameters  as  angles  of  arrival,  number  of  signals,  etc. 
Schmidt  [28]  and,  independently,  Bienvenu  [30]  were  the  first  to  correctly 
exploit  the  measurement  model  in  the  case  of  a  sensor  array  of  arbitrary 
shapes.  Schmidt's  algorithm,  called  MUSIC  (Multiple  Signal  Classification), 
identifies  two  distinct  eigenspaces.  The  space  associated  with  the  smallest 
eigenvalue  which  appears  with  a  multiplicity  of  (m-d),  where  m  is  the  num¬ 
ber  of  sensors  and  d  is  the  number  of  targets,  assuming  m>d,  is  called  the 
noise  subspace.  The  space  associated  with  the  d  non  zero  eigenvalues  is 
called  the  signal  subspace.  The  angles  of  arrival  of  the  sources  are 
estimated  via  a  search  procedure  which  consists  of  choosing  directional 
vectors  and  correlating  them  with  the  noise  space  generated  by  the  noise 
eigenvectors  corresponding  to  the  noise  eigenvalues.  Since  the  noise  space 
is  orthogonal  to  the  signal  space,  the  angles  of  arrival  are  the  peaks  in 
the  reciprocal  of  this  correlation.  Computationally,  this  algorithm  is  very 
inefficient.  Another  disadvantage  of  MUSIC  is  that  it  cannot  handle  com¬ 
pletely  correlated  sources.  A  pre-processing  technique  called  Spatial  Smoo¬ 
thing  [31],  was  then  suggested  for  the  case  of  a  linear  uniformly  spaced 
array.  Spatial  Smoothing  uses  a  set  of  L  contiguous  subarrays  (L<m). 

(m-L+1)  correlation  matrices  are  added  up  to  form  a  new  correlation  matrix. 
This  matrix  is  shown  to  be  non  singular.  However,  this  processing  decreases 
the  effective  aperture  size  which  reduces  the  ability  of  the  array  to 
detect  a  sizable  number  of  sources  since  the  two  are  directly  related. 

Non  search  procedures  were  introduced  to  reduce  such  problems  as 
computational  complexity,  storage,  etc.,  which  are  inherent  in  search  pro¬ 
cedures.  A  non  search  procedure  using  MUSIC  was  initiated  by  Barabell  [32] 
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Cor  the  case  of  a  linear  uniformly  spaced  array.  Instead  of  plotting  the 
MUSIC  spectrum,  a  root  finding  procedure  is  developed,  where  the  roots  of 
the  polynomial  are  functions  of  the  parameters  of  interest.  An  improvement 
has  been  noticed,  especially  for  the  case  of  closely  spaced  targets. 

Pisarenko's  method  [33]  is  based  upon  the  fact  that  tne  covariance 
matrix  has  a  smallest  eigenvalue  of  multiplicity  (m-d).  Thus,  a  repeated 
eigenanalysis  is  required  to  determine  the  multiplicity  of  the  smallest 
eigenvalue.  In  practice,  however,  because  the  correlation  matrix  is 
evaluated  from  finite  data  samples,  it  is  very  difficult  to  count  the  mul¬ 
tiplicity  of  the  smallest  eigenvalue.  More  sophisticated  approaches  based 
on  statistical  considerations  have  been  developed  [34-38].  Yet,  the  choice 
of  a  subjective  threshold  make  these  methods  undesirable.  Wax  and  Kailath 
[39]  developed  methods  based  on  the  information  theoretic  criteria  intro¬ 
duced  by  Akaike  (AIC)  and  by  Schwartz  and  Rissanen  (MDL) .  The  number  of 
signals  is  derived  by  minimizing  the  MDL  or  AIC  functions.  Recently,  Zhao 
et.  al.  [40]  showed  that  AIC  is  inconsistent  and  further  suggested  a  family 
of  consistent  estimators  of  which  MDL  is  a  member. 

ESPRIT  (  Estimation  of  Signal  Parameters  via  a  Rotational 
Invariance  Technique)  was  later  proposed  by  Roy  and  Kailath  [43].  It  is 
based  on  an  invariance  introduced  by  the  array  geometry.  It  is  shown  to  be 
robust  and  computationally  very  efficient.  The  angles  of  arrival  of  the 
sources  are  directly  related  to  the  generalized  eigenvalues  of  a  matrix 
pencil  formed  from  the  data.  However,  the  estimates  obtained  v* th  this  al¬ 
gorithm  are  very  biased  at  low  SNR.  Although  the  original  version  of  ESPRIT 
employed  a  least  squares  criterion,  a  total  least  squares  procedure  [44,45] 
was  then  devised  to  overcome  this  disadvantage. 

Yet,  ESPRIT  can  not  be  used  in  the  case  of  completely  correlated 
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sources.  H.  Ouibrahim  has  shown  that  ESPRIT  is  only  one  of  several  possible 
operators  that  could  be  used  in  the  Matrix  Pencil  Approach  [1].  Two  other 
operators  have  been  proposed  which  are  referred  to  as  the  Summation  Opera¬ 
tor  and  Moving  Window  operator.  The  Moving  Window  operator  is  shown  to  hold 
even  in  the  case  of  fully  correlated  sources,  thus  outperforming  ESPRIT. 

Ouibrahim  also  showed  [46 J  that  the  Moving  Window,  Prony's  method 
and  Pisarenko's  method  are  all  equivalent  in  the  sense  that  they  check  the 
dependence/ independence  of  some  data  vectors.  S.  Mayrargue  [47,48]  showed 
that  ESPRIT,  TAM  (Toeplitz  Approximation  Method)  [49,50]  and  Tufts  and 
Kumaresan's  method  [51,52]  are  all  equivalent  in  the  sense  that  they  all 
solve  the  same  multidimensional  system. 

All  of  this  work  was  done  for  the  case  of  narrowband  signals.  The 
case  of  wideband  signals  is  much  more  difficult  and  has  been  studied  by 
only  a  few  people.  Su  and  Morf. [53,54]  suggested  using  a  Modal  Decomposi¬ 
tion  of  the  signals  along  with  MUSIC  to  solve  for  the  angles  of  arrival.  A 
set  of  angles  of  arrival  is  obtained  for  each  pole  in  the  received  spec¬ 
trum.  Thus,  more  signal  parameters  like  poles  and  residues  have  to  be 
estimated. 

A  different  way  of  approaching  the  problem  is  to  decompose  the 
wideband  signal  into  narrowband  signals  and  then  use  the  well  known  narrow- 
band  algorithms.  At  this  stage  two  schemes  were  developed;  post  averaging 
schemes  and  pre-averaging  schemes. 

The  first  scheme  was  used  by  Wax  et.  al.  i 55 ] .  They  suggested  a 
modified  version  of  the  MUSIC  algorithm.  The  spatial  spectral  estimate  was 
formed  by  averaging  (either  geometrically  or  arithmetically)  the  MUSIC  in¬ 
ner  products  for  each  of  the  frequencies  considered.  This  method  is  termed 
as  Incoherent  processing.  It  cannot  be  employed  in  the  case  of  signals  ex- 
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periencing  multipath. 

Vang  and  Kaveh  [56]  proposed  a  pre-processing  scheme  using  linear 
transformations  to  combine  the  various  sub-bands  into  one  single  band.  This 
is  called  Coherent  Wide-Band  (CVB)  processing.  Only  one  eigendecomposi tion 
is  then  performed  at  this  reference  frequency.  This  algorithm  outperforms 
the  previous  one  in  many  ways  such  as  computation,  applicability  in  the 
case  of  coherent  sources,  etc.  However,  it  suffers  from  the  fact  that 
preliminary  estimates  of  the  angles  of  arrival  are  needed  in  order  to  form 
the  linear  transformations.  If  these  angles  are  clustered  within  a  beam- 
width,  then  the  method  performs  well.  Otherwise,  spatial  prefiltering  is 
needed.  . 

Another  wideband  method  was  introduced  by  Buckley  and  Griffiths 
[58]  called  BASS-ALE  (Broad-Band  Signal  Subspace  Spatial  Spectrum  Estima¬ 
tion).  They  generate  a  signal  subspace  using  a  "stacked",  vector  snapshot. 
The  vectors  obtained  at  each  subband  are  stacked  on  top  of  one  another. 

This  algorithm  is  computationally  more  expensive  than  CVB. 

A  new  coherent  wide  band  algorithm  was  proposed  by  Kumaresan  and 
Shaw  [60].  They  make  use  of  a  bilinear  transformation  to  transform  each 
sub-band  into  the  reference  band.  This  can  be  considered  as  a  one  step  al¬ 
gorithm.  Also  it  does  not  need  the  a  priori  knowledge  of  the  angles  of  ar¬ 
rival.  However,  it  can  o.:ly  be  applied  to  electronically  small  arrays.  Its 
performance  deteriorates  for  targets  near  the  endfires. 

Recently,  Ottersten  and  Kailath  [61,62]  proposed  extending  the 
ESPRIT  algorithm  to  wide  band  signals.  They  use  the  same  model  as  is  de¬ 
scribed  in  Su  and  Morf  [36]  and  ESPRIT  is  applied  to  determine  the  source 
locations  for  each  of  the  signal  poles. 
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1.4  WORK 


OUTLINE 


The  remainder  of  this  work  is  organized  as  follows.  Chap¬ 

ter  2  introduces  the  principles  of  signal  subspace  processing,  the  pencil 
theorem,  (which  is  the  basis  of  this  work),  and  the  Matrix  Pencil  approach. 
New  ways  for  evaluating  the  generalized  eigenvalues  are  also  proposed.  All 
of  the  algorithms  discussed  above  assume  an  ideal  sensor  environment  in 
which  each  sensor  is  assumed  to  exist  by  itself.  Therefore,  effects  of 
reradiation  from  these  sensors  are  completely  neglected.  Chapter  3  deals 
with  the  problem  of  mutual  coupling  between  the  sensors.  Compensation  tech¬ 
niques  are  developed  for  the  Matrix  Pencil  Approach  using  the  Moving  Window 
and  the  ESPRIT  operators.  Chapter  4  is  devoted  to  the  extension  of  the 
Moving  Window  to  the  broadband  case.  Three  methods  are  devised.  The  first 
is  original  in  the  sense  that  the  Matrix  Pencil  approach  is  utilized  with  a 
signal  model  not  used  previously  in  other  approaches.  The  signals  are 
identified  by  their  poles  (natural  frequencies)  and  their  angles  of  ar¬ 
rival.  The  second  method  utilizes  the  same  model  used  by  Su  and  Morf.  How¬ 
ever  the  analysis  is  carried  out  completely  in  the  time  domain.  The  third 
makes  use  of  the  CWB  of  Wang  and  Kaveh.  Chapter  5  analyzes  the  effects  of 
the  noise  and  perturbations  due  to  sensor  spacing.  A  measure  is  introduced 
and  a  geometric  upper  bound  is  derived  for  the  Moving  Window  and  ESPRIT  op¬ 
erators.  Conclusions  and  recommendations  for  future  research  are  presented 
in  chapter  6. 
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CHAPTER  2 


SIGNAL  SUBSPACE  PROCESSING 


2.1  PRINCIPLE 

Assume  there  are  d  sources  emitting  signals  s^t)  ;k=l,2,.  .  .  ,d, 
which  are  impinging  on  a  linear  array  composed  of  m  sensors.  It  is  assumed 
that  d<m.  Let  the  superscript  T  denote  transpose.  The  received  signal  at 
the  i-th  sensor  can  be  modeled  as 

d 

Xi(t)  =  I  slc(t)ai(ek)  +  n^t),  (2.1-1) 

k-1 


where 


ai(®k)  *s  t*ie  relative  response  of  the  i-th  sensor  to  the  k-th 

source, 

Sk(t)  is  the  complex  envelope  of  the  signal  emitted  by  the  k-th 


source, 

n^  is  the  additive  noise  assumed  to  have  zero  mean  and  variance  . 
In  vector  notation,  this  can  be  written  as 

X  *  A  S  +  N  (2.1-2) 


where 


XT  ■  [x^,X2»  •  •  •  ,xm]-  (mxl)  received  signal  vector, 

ST  »  [s^,S2»  •  •  •  ,sd]«  (dxl)  impinging  signal  vector, 

Nt  «  (n^,n2»  •  •  •  ,nm]»  (mxl)  noise  vector, 

A  -  [a^  §2  •  •  •  ad  1*  (mxd)  direction  matrix, 

a^  »(a^(0j)  a2(0i)  .  .  ,am^®i^J*  (mxl)  i£^  direction  column 

vector  of  A. 

In  all  the  subspace  approaches  that  have  been  proposed  the  signals 
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and  noise  are  assumed  to  be  statistically  independent  and  the  noises  nj  are 
assumed  to  be  independent  from  sensor  to  sensor  with  a  correlation  matrix 
given  by  a2!  where  I  is  the  identity  matrix.  Let  the  superscript  H  denote 
the  Hermitian  Transpose.  The  spatial  covariance  matrix  is 
R  -  E[X  XH]  =«  E[(AS+N)(AS+N)h] 

-  E[ AS  SflAH]  +  E{N  Nh} 

-  AE[S  Sh]Ah  +  a2!  (2.1-3) 

Let  S«E{S  S«J.  Then  R  can  be  written  as 

R  =.  ASAH  +  a2!  (2.1-4) 


where 


R  is  an  (mxm)  matrix. 

Let  >  A2  >  A3  >,  .  .  •  ,>>„,}  be  the  set  of  eigenvalues  of  R.  Let 
{Zi*!l2*¥3’  *  *  *  * Ym  )  be  the  set  o£  the  corresponding  eigenvectors. 

If  S  is  non  singular  and  with  the  assumption  that  m  >  d,we  can 


show  that 


1)  the  minimum  eigenvalue  of  R  is  g2  with  multiplicity  (m-d),i.e, 

xd+l*xd+2"xd+3**  *  *  '  =*Xm”Xmin=<T2* 

2) the  eigenvectors  associated  with  the  minimum  eigenvalue,  V<j+i> 

¥d+2»  ^d+3*  •  •  •  are  orthogonal  to  the  space  spanned  by  the  columns 

of  A. 

These  results  lead  to  the  following  direction  finding  approach  : 

1)  determine  the  number  of  sources  d  from  the  multiplicity  of 


\iin* 

2)  use  the  orthogonality  relation  between  the  direction  vectors  of 
the  impinging  sources  and  the  eigenvectors  corresponding  to  Amin  to  yield 
the  directions  of  arrival  of  the  sources.  We  just  have  to  "search"  for 
those  direction  vectors  that  are  orthogonal  to  the  eigenvectors  cor- 
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responding  to  X^j^.  For  this  reason  methods  based  on  this  approach  are 
called  search  procedures. 

2.2  PENCIL  THEOREM 
2.2.1  Theorem 

Denote  by  C  the  field  of  all  complex  numbers. Consider  two  matrices 
M  and  N  of  size  (kxp).The  set 

{  M-XN  {  X  e  C  } 

is  said  to  be  a  matrix  pencil.  The  matrices  M  and  N  are  required  to  have 
the  following  decompositions 
M  «  E  F 
N.EDF 

where 

E  is  a  (kxd)  matrix  and  k  >  d 
F  is  a  (dxp)  matrix  and  p  >  d 

D  is  a  (dxd)  diagonal  matrix  where  d^  denotes  the  i-th 
diagonal  element. 

If  M  and  N  are  two  matrices  which  have  the  decompositions  cited 
above  and  if  E,  F  and  D  are  all  of  rank  d,  then  the  rank  of  the  matrix  pen¬ 
cil  M-XN  is  decreased  by  1  whenever 

Xi  *  (dii)  ^  ;  i*l»2,...fd. 

These  values  of  Xj  are  known  as  the  generalized  eigenvalues  of  the  matrix 
pencil. 

Proof 

Since 

M-EF 

and 
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Thus, 


N-EDF, 

M-XN  -  EF-XEDF  -  E(I-XD)F. 


rank  (M-XN)  =  rank(E(I-XD)F)  =*  min{rank(E) , rank(F) , rank(I-XD)} . 
However,  by  assumption 

rank(E)»rank(F)*d 

and 

rank(I-XD)  is  of  rank  d  as  long  as  (1-X}d^)*0  ;i=l,2,.  .  ,d. 

If  (1-Xjdjj)*0,  which  implies  that  Xj=(djj)-1  , rank(I-XD)=d-l .  Therefore, 
the  rank  of  (M-XN)  is  reduced  by  1  whenever 

.  Xj=(dj^)  ^  j  i*l , 2 d . 


2.2.2  Determination  of  rank  reducing  values 

Note  from  above  that  two  cases  may  occur.  If  k=p,  M  and  N  are 
square  matrices.  The  set  of  the  generalized  eigenvalues  of  the  pencil  M-XN 
is  defined  to  be  the  set  of  all  elements  Xj  such  that 

det(M-XiN)=0. 

When  the  generalized  eigenvalues  are  distinct, the  rank  of  M-XN  is  reduced 
by  1  whenever  X  equals  one  of  these  values.  In  the  case  where  k>p,  M  and  N 
are  non  square  matrices. Det(M-XjN)  no  longer  exists  since  the  pencil  is  not 
square.  For  this  reason  we  have  to^make"  the  pencil  matrix  a  square  one. 
This  can  be  done  by  premultiplying  the  pencil  M-XN  by  either  MH  or  NH.  This 
can  also  be  done  by  postmultiplying  the  pencil  M-XN  by  either  MH  or  N^.  We 
obtain 

MH(M-XN)  -m=m-xm% 
or 

Nh(H-XN)  =NhM-XNhN. 
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MH(M-XN)  and  NH(M-XN)  are  both  square  matrices  (pxp).  Notice  that 
MH(M-XN)«(EF)H(EF-XEDF)»PHEHEF-XFHEHEDF 
=.FHEhE(I-XD)F 

and 

NH ( M- XN ) . ( EDF ) H ( EF- XEDF ) = FHDHEHEF- XFHDHEHEDF 
*FHDHEHE(I-XD)F. 

Both  equations  have  the  decompositions  required  by  the  pencil  theorem 


since 


F^E^E  and  F  are  of  rank  d 


FHDHEHE  and  F  are  of  rank  d. 

Because  (I-XD)  arises  in  all  of  these  equations,  ve  can  say  that  the  rank 
reducing  values  of  the  pencils  Mh(M-XN)  and  N^(M-XN)  are  identical  to  those 
of  the  pencil  (M-XN). 

Note  that  the  pencils  MH(M-XN)  and  NH(M-XN)  have  p  generalized 
eigenvalues.  However,  there  is  a  method  which  relies  on  a  singular  value 
decomposition  (SVD)  of  the  matrices  M  and  N  to  obtain  a  matrix  pencil  which 
has  exactly  d  generalized  eigenvalues  which  in  turn  are  the  rank  reducing 
values  of  the  pencil  (M-XN).  The  singular  value  decompositions  of  the 
matrices  M  and  N  result  in 


where 


M-lWm* 

N-U„S„VnH, 


Um  .  [Uni!  Umj 


Umj  »  i-th  right  singular  vector  of  M  =  (kxl), 


(2. 2. 2-1) 


Urm.  ]  «  (kxk) 


I  0 


(kxp) 
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Eju  *  diag  {  cm j  om2  •  •  •  )  *  (dxd), 

onij  =  i-th  singular  value  of  M 
vm  *  IY”1  Vra2  .  •  •  Vrap  ]  «  (pxp) 

Vmj  =  i-ch  left  singular  vector  of  M  =  (pxl), 
Un  =  [Unj  Un?  .  .  .  Uni,  ]  =  (kxk) 

Unj  =  i-th  right  singular  vector  of  N  =  (kxl), 


S 


n 


'  I  0  ' 

.  0  |  0  . 


(kxp) 


En  =  diag  {  cm ^  on2  •  •  •  cmd  }  =  (dxd), 
on^  *  i-th  singular  value  of  N 

vn  *  fYSl  ¥H2  •  •  ♦  ^*lp  1  *  (PXP> 

Vni  *  i-th  left  singular  vector  of  N  *  (pxl). 

Collecting  the  d  principal  left  and  right  singular  eigenvectors  (cor¬ 
responding  to  the  d  non  zero  singular  values),  we  have 
"  -  M*  >  Um'  H, 

and  (2. 2. 2-2) 

N  -  N«  -  U/  Zn  (Vn*)H, 

where  t  denotes  truncated  and 

Um1  -  (U51  Um2  .  .  .  Um^  J  *  (kxd) 

Vm 1  -  [Vni  Vm2  .  .  .  Vmd  ]  -  (dxp) 

un*  *  [HUl  MH2  •  *  *  Mild  J  *  (kxd) 

Vnl  “  [Vn^  Vn2  .  .  .  Vnd  ]  *  (dxp). 

Therefore, 

(M-XN)  »  (M^-XN* )  -  (UmtEm(Vmt)H-XUnt!:n(Vnt)H).  (2. 2. 2-3) 
Pre-multiplying  and  post-multiplying  both  sides  of  equation  (2. 2. 2-3)  by 
(U,,*)*1  and  Vnl  respectively,  we  get 
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(Uj,*)3  (Mt-XN^Vn1  =  ((Unt)HUmtEm(Vmt)%nt  -**n>-  (2. 2. 2-4) 

This  pencil  is  square  and  is  of  dimensions  (dxd)  and  thus  has  d  rank  reduc¬ 


ing  values  which  are  exactly  its  generalized  eigenvalues. 


2.3  MATRIX  PENCIL  APPROACH 

As  was  stated  earlier,  the  matrix  pencil  approach  is  based  upon 
the  pencil  theorem.  Two  matrices  are  formed  from  the  data  generated  at  the 
outputs  of  the  sensor  array.  The  generalized  eigenvalues  of  the  pencil  thus 
formed  are  shown  to  contain  the  information  about  the  locations  of  the 
targets.  In  this  section  we  discuss  two  different  operators. 

2.3.1  ESPRIT 

Consider  a  planar  array  of  arbitrary  geometry  composed  of  2m 
sensors  arranged  in  pairs  so  as  to  form  m  doublets  having  the  same  direc¬ 
tional  orientation  with  respect  to  each  other.  The  elements  of  each  doublet 
have  identical  directional  g,.in  patterns  and  are  translationally  separated 
by  a  known  displacement  A  (Fig.  2-1).  Other  than  the  obvious  requirement 
that  each  sensor  have  non-zero  gain  in  the  directions  of  the  emitting  sig¬ 
nals,  the  beam  pattern  of  the  elements  in  the  doublet  are  totally  ar¬ 
bitrary.  Assume  there  are  d  (m>2d)  narrowband  sources  centered  at  frequency 
«,  and  that  the  sources  are  located  sufficiently  far  from  the  array  such 
that  the  wavefronts  impinging  on  the  array  are  planar.  Assume  the  sources 
are  located  at  azimuthal  angles  9^,  k-1,  2,  .  .  .,d  and  emitting  signals 
whose  complex  envelopes  are  denoted  by  s^,  k»l,  2,  .  .  .,  d.  Additive  noise 
is  present  at  all  2m  sensors  and  is  assumed  to  be  stationary  with  zero- 
mean  and  variance  a2.  The  signals  received  at  the  two  sensors  in  the  i-th 
doublet  can  be  expressed  as 
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(2.3-1) 


d 

vj ( t )  -  Z  s^Og^ 9k)+nix(t)  =  xi(t)+nix(t) 
k»l 


d 

v^t)  -  Z  gi(©it>+niy^t>  ■  yi(t)+niy(t)f 

k-1 


where  g^G^)  is  the  gain  response  of  the  i-th  sensor  to  a  source  arriving 
at  angle  9^. 

Two  vectors  V  and  V  are  then  formed  where 
VT  -  [Vl  v2  v3  .  .  .  vm] 

and 

WT  =  [wx  w2  w3  .  .  .  wm]. 

V  and  tf  can  be  written  as 


where 


V  =  GS  +  N. 


X  +  N* 


(2.3-2) 


V  a  GSS  +  Ny  *  Y  +  Ny  , 

G,  *,  S,  Njj  and  Ny  are  the  following  matrices: 

G  “  Ifil  82  '  *  •  Sd  J=  (mxd)  gain  matrix, 

gi  *  [«i( ©i >  ®2< ®i )  •  •  •  gm(9i>]=*(n«l)  ith  gain  column  vector  of  G. 
ST  »  (s^  s2  •  •  •  s<j]»  (dxl)  impinging  signal  vector, 

5xT  *  lnlx  n2x  •  *  •  nraxl*  (®xl)  noise  vector, 

NyT  ■  [n^y  n2y  .  .  .  nmy]-  (mxl)  noise  vector, 


and 


*  -  diag{eJ*l  eJ*2  .  .  .  e^d) 

where 


♦k*-(«A/c)sin(ek)  ,  k-1,  2,  .  .  .,  d.  (2.3-3) 

Assuming  that  the  signals  and  noise  are  statistically  independent  and  that 
the  noise  components  are  uncorrelated  from  sensor  to  sensor,  we  have 
E[V  VH]  -  GSGh  +  a2  I 
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(2.3-4) 


E[ V  WH1  -  GS^G11 

where  I  is  the  (mxm)  identity  matrix  and  S*E[S  SH].  Consider  the  matrices  M 
and  N,  where 

M  -  E[ V  VH]  -  o2!  -  E[X  XH]  -  GSGH 

(2.3-5) 

N  -  E[V  VH]  =.  E[ X  YH]  =  GS^G0  . 

Consider  the  pencil  M-XN. 

M-XN  -  (GSGH)-X(GS*HGH)  *  GS(I-X*H)GH.  (2.3-6) 

In  the  case  where  the  sources  are  not  fully  correlated  (hence,  S  is  not 
singular)  and  the  direction  of  arrival  of  the  sources  are  all  distinct,  M 
and  N  are  of  rank  d.  The  rank  of  this  pencil  is  reduced  by  1  whenever 

Xj^eJ^k  *  exp{-j(wA/c)sin(0jt)}  ;  k»l,2,.  .  .,d.  (2.3-/) 

The  angles  of  arrival  are  then  given  by 

9k.»sin-1{jcln(Xk) /(«£)}  ;  k=l,2,  .  .  .,d.  (2.3-8) 

2.3.2  MOVING  WINDOW 

Assume  we  have  a  linear  array  composed  of  m  identical  sensors  with 
uniform  spacing  D.  Let  there  be  d<  m  narrowband  sources  located  at 
azimuthal  angles  0^  ;  k*l,  2  ,  .  .  d  (Fig.  2-2).  Let  the  complex  en¬ 
velopes  of  the  emitting  signals  be  denoted  by  sjc(  t ) ;  k»l,  2,  .  .  .,  d.  As¬ 
sume  the  sources  are  in  the  far  field  such  that  planar  wavefronts  arrive  at 
the  array.  With  reference  to  the  first  sensor,  the  received  signal  at  the 
ith  sensor  is  modeled  as 

d 

vi(t,9)«E  Sk(t)ai(©k)+ni(t):*xi+ni;  i*l»  2,.  .  .  ,m  (2. 3. 2-1) 

k-1 

where 

ai(©k)  =  a^  e^i_1^^k 
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=  -wDsinCO^J/c  ;  k=l,  2,  .  .  .  ,d, 

(O  is  the  center  frequency  of  the  plane  waves 
c  is  the  propagation  speed  of  the  waves 
0  is  the  sensor  spacing 

ajcaa(6lc)  is  the  beam  pattern  in  the  direction  of  the  k-th  emitter 

and 

n^(t)  is  the  additive  noise  assumed  to  be  zero-mean. 

2.3.2. 1  Original  Version 

in  the  original  formulation  of  the  matrix  pencil  approach  using 
the  moving  window,  Ouibrahim  [1]  generated  (d+1)  vectors  of  length  (m-d). 
These  vectors  were  then  used  to  form  two  (dxd)  matrices  M  and  N.  Unlike  the 
ESPRIT  approach,  these  matrices  are  of  rank  d  even  in  the  presence  of  fully 
correlated  sources.  The  only  restriction  is  that  all  sources  have  distinct 
angles  of  arrival.  For  this  case,  there  will  always  be  d  generalized  eigen¬ 
values  of  the  pencil  (M-XN)  regardless  of  the  size  of  the  sensor  array.  In 
the  absence  of  noise,  Ouibrahim's  choice  results  in  generalized  eigenvalues 
whose  corresponding  eigenvectors  span  the  signal  subspace.  However,  in  the 
presence  of  noise,  the  generalized  eigenvalues  are  contaminated  by  the 
noise  so  that  neither  the  signal  nor  the  noise  subspaces  can  be  identified. 
This  explains  why  the  method  performs  badly  t jr  signal  to  noise  ratios 
smaller  than  15  dB.  To  overcome  this  deficiency,  it  is  preferable  to  choose 
a  window  of  length  L<(m-d)  and  then  form  (m-L+1)  vectors  Vn  of  length  L 
where 

ZnT  *  lvn>  vn+l  *  •  •  •»  vn+L-lI  5  n=1>-  •  • , (m-L+1). 

The  limits  of  L  will  be  derived  later.  It  can  be  shown  that  Vn  can  be  writ¬ 
ten  as 
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(2. 3. 2. 1-1) 


Yn  -  +  N„  -  Xn  +  Nj 

where 


■  1  1. 

ej  +1  e  j  *2 

•  •  • 

.  1. 

.  ej  4*d 

•  •  • 

_  ej(L-i)+i  ej(L-D+2 ; 

;  ;j(L-D#d 

#  =>  diag{ej^l  ej^2  .  .  .  e^d} 

B  »  diag{alf  a2,  •  •  • .  ad  } 

S  «  t S1 ’ s2 )* 

m 

Nn  ■  Inn>nn+1» *  *  *  * »nn+L-l  1* 

The  two  matrices  and  are  then  formed 


- 

- 

_ i 

'  V1  •  •  •  V(m-L) 

1  1 

v2  *  •  *  V(m-L+1) 

Mi  m 

Yl  *  *  *  Y(m-L) 

s 

•  • 

•  • 

1  1 

.  4  4  J 

•  • 

■  VL  •  •  •  v(m-l) 

'  T  T 

'  v2  •  •  •  V(m-L+1) 

Ni  - 

1  1 

Y2  •  •  •  Y(m-L+1) 

v3  •  •  •  v(m-L+2) 

•  • 

•  • 

•  • 

1  1 

.4  4 

•  • 

-  V(L+1)  •  •  •  vm 

Using  the  expression  for  Yn  I  n-l»  2 . (m  -L+l),  the  matrices  Mj 

become 

Mi  -  [  AiBS  A!B*S  .  .  .  A1B*<ra-L-1>S  ]+  [  Nj  N2  .  .  .  N(m_L)  ], 
NX  -  [  AiBfS  A1B*2S.  .  .  A1B*<ra-L-1>S  J+  [  N2  N3  .  .  .  N(m_L+1) 
This  can  also  be  written  as 

MX  -  AXB  [S  #S  .  .  .  *("J-L-l)s  ]+  [  Ni  N2  .  .  .  N(m_L)  ], 

NX  -  AiBt  [S  *S  .  .  .  *<m-L-l)s  ]+  [  N2  N3  .  .  .  N(m_L+1)  J. 


3.2. 1-2) 


3.2. 1-3) 


and  Nj 


]• 
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Let  C  be  the  matrix 


C«[  S  *S  .  .  .  *<m-L>s  J. 

It  can  be  shown  that  C  can  be  expressed  as 
C  -  D  U, 

where 

D*diag{S!  s2  .  .  .  sj} f 

1  ej*l  .  .  .  ej<™-L-l)+i 

1  ej*2  .  .  .  ei<m-L-1)*2 

•  •  • 

U  -  .  . 

.  i  ij+d  .  .  .  ej<m-L-1)+d 

Let  N'  and  Nn  be  the  matrices 

N'  -  (  NX  N2  .  .  .  N(m_L)  ], 

N"  -  [  N2  N3  .  .  .  N(m_L+1)  ]. 

Then 

Mx  -  A1BDU+N' 

and  (2. 3. 2. 1-4) 

Nx  -  A1BD«J+N". 

Assuming  that  the  signals  and  noise  are  statistically  independent  and  that 
the  noise  components  are  uncorrelated  from  sensor  to  sensor  with  zero  mean 
and  variance  a2,  we  get 

ElM^l-^VU  +  L  a2  I(m_L)  (2. 3.2. 1-5) 

EIN^M^-U8*0™  +  L  a2  Ii(m_L)  (2. 3. 2. 1-6) 

where  I(m_L)  is  the  (m-L)x(m-L)  identity  matrix,  Ii(m_L)  is  the  matrix 

'0100.  .  .  .0 
0  C  1  0  ....  0 

0  0  0  1  ....  0 

•  •  •  •  • 

Il(m-L)  "  •  •  *  * 

•  •  •  •  • 

0  0  0  0  ....  1 

.0000.  .  .  .0 
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and  V  is  the  matrix  V»E[DHBHA1HA^BD] .  Defining 

Ppq-t  ^(i-D(W. 

i*l 

Spq  *  Ft  sqsp  ]  i 
apq  =  aJsp, 

the  matrix  V  can  be  written  as 

'  sllallFll  s21a21F21  •  •  •  •  sdladlFdl 

s12a12F12  s22a22F22  ••  •  sd2ad2Fd2 

•  •  • 

V  - 

•  •  •  • 

■  sldaldFld  s2da2dF2d  •  •  •  •  SddaddFdd  . 

Define  the  matrices  M  and  N  such  that 

M  -  EIM^]-  L  a2  I(m_L)  =UHVU  (2.3.2. 1-7) 

N  -  L  a2  Ii<m_L)»UH«®VU.  (2.3.2. 1-8) 

Note  that  M  and  N  are  (m-L)x(m-L)  matrices.  Assuming  that  (m-L)  >  d,  M  and 
N  will  always  be  of  rank  d.  Therefore,  the  limits  on  L  are 

d  <  L  <(m-d) . 

Consider  the  matrix  pencil  M-AN. 

M-AN  -  UHVU-AUH*HVU  -  UH(I-X*H)VU  (2.3.2. 1-9) 

which  satisfies  the  requirements  of  the  pencil  theorem.  Hence,  the  values 
of  A  for  which  the  rank  of  M-AN  decreases  by  1  are  given  by 

\  m  e-3*k  ;  k*l,2,...,d.  (2.3.2.1-10) 

The  angles  of  arrival  are  given  by 

\  -  sin~1{jcln(Alc)/wD} ;  k«l,2,...,d.  (2.3.2.1-11) 

2. 3. 2. 2  Nev  Version 

From  equations  (2. 3. 2. 1-2)  and  (2. 3. 2. 1-3),  we  see  that  we  can 
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form  (L+l)  vectors  Zr  of  length  (m-L)  where  Zr  is  given  by 
ZrT  m  {Vj.  ,  Vj.+j^  »••••»  Vp+m_L_i}  *  r-1, .  .  •  ,(L+1), 
and  *  denotes  complex  conjugate.  We  have  purposely  chosen  to  take  the  con¬ 
jugate  so  that  we  deal  with  correlation  matrices  of  the  form  E[Z  Z^] .  Note 
that  this  is  the  usual  form  for  the  correlation  matrix  of  the  vector  Z.  it 
follows  that  M}  and  Nj  are  simply 


r  -  ZiH  i 

-  I2H  - 
• 

r  <-  z2h  -»  i 

-I3H  - 

* 

and  Nj  - 

.  -  k"  -  ■ 

Z(L+1)H  -*  . 

Therefore, 


E  [M^] 

and 

EIN^] 


L 

£ 

i-1 

L 

£ 

i-1 


E[ZiZiHJ 

EIZ(i+l)ZiH]. 


(2. 3. 2. 2-1) 


(2. 3. 2. 2-2) 


The  advantage  of  formulating  the  matrix  pencil  approach  using  the  moving 
window  in  this  fashion  will  be  seen  later  when  mutual  coupling  is  present 
at  the  sensor  array.  Using  the  original  version,  it  was  necessary  to  employ 
a  minimum  mean  squared  error  estimation  scheme  to  compensate  for  the 
mutuals.  Using  the  new  version,  it  was  possible  to  devise  a  scheme  that 
provided  a  direct  compensation  for  the  mutuals. 


3.  Computer  Simulation  of  the  Moving  Yindov  Approach 

The  model  used  in  the  computer  simulation  consisted  of  two  in¬ 
coherent  sources  (d-2)  incident  on  a  linear  array  of  eight  uniformly  spaced 
sensors  (m-8).  For  convenience,  the  sensors  are  assumed  to  be  omnidirec¬ 
tional.  The  Rayleigh  resolution  of  this  array  is  given  by 
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2/(m-l) (180/re)  =16.37°. 

The  2  sources  were  assumed  to  be  located  at  dj.*16°  and  @2»24°.  The  angular 
separation  o£  these  sources  is  8°  which  is  half  the  Raleigh  resolution.  The 
amplitudes  of  the  sources  were  generated  as  statistically  independent  com¬ 
plex  random  variables  with  zero  mean  and  unit  variance.  The  noise  was  also 
chosen  to  be  complex  Gaussian  with  zero  mean  and  unit  variance.  The  sensors 
were  positioned  at  half  wavelength  so  that  (wgD/c)=Tt.  100  snapshots  were 
taken  for  each  of  the  50  runs.  The  length  of  the  window  was  varied  frem 
L-d-2  to  L*(m-d)»6.  In  the  plots  of  the  sample  variance  and  the  mean- 
squared  error,  the  y-axis  is  defined  as 
y*10  log10(-). 

Let  9^  be  an  estimate  of  9  obtained  at  the  k-th  run  (K  is  the  number  of 
runs).  The  sample  mean  (ME),  the  sample  variance  (Var)  and  the  mean-squared 
error  (MSE)  are  defined  respectively  as 

K 

ME(9)  .  (1/K)  Z  9k, 
k-1 

K 

Var (9)  -  (1/K)  Z  (9k-ME(9))2, 
k-1 

K 

MSE(9)  .  (1/K)  Z  (9k-9)2. 
k-1 

The  results  are  shown  in  Fig.  2.3  to  2.8.  From  figures  2.3  and  2.4,  note 
that  above  15  dB,  the  choice  of  L  is  not  important  as  far  as  the  sample 
mean  is  concerned.  However,  below  15  dB,  the  performance  of  the  Moving 
Vindow  degrades  considerably  with  the  choice  of  L-5  and  L=6  (  Ouibrahim's 
choice).  For  L-2,  3  or  4,  the  performance  of  the  method  is  comparable.  By 
choosing  a  windov  of  length  L-2  (smallest  possible),  an  improvement  of  al- 
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most  23  dB  is  achieved  at  a  signal  to  noise  ratio  of  5  dB  in  the  sample 
variance  and  the  mean-squared  error.  This  improvement  is  due  mainly  to  the 
recognition  of  both  subspaces  (signal  and  noise)  and  the  use  of  the 
singular  value  decomposition  (SVD).  The  SVD  is  known  to  be  robust  even  in 
the  case  of  very  ill  conditioned  matrices.  The  reduction  of  the  2  matrices 
involved  in  the  matrix  pencil  to  the  desired  dimension  (d)  through  the  SVD 
allowed  us  to  effectively  use  the  IMSL  routine  (EIGZC)  which  computes  the 
generalized  eigenvalues. 
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SAMPLE  MEAN 


ANGLE=16° 


Fig.  2.3  Sample  Mean  of  the  Angle  Estimate 
at  16° 
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SNR  (dB) 


Fig.  2.6  Sample  Variance  of  the  Angle  Estimate 
at  24°. 
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MEAN-SQUARED  ERROR  (dB) 


MEAN-SQUARED  ERROR  (dB) 


CHAPTER  3 


COMPENSATION  FOR  MUTUAL  COUPLING 

In  the  methods  considered  earlier,  each  sensor  was  treated  as  if 
it  existed  by  itself.  However,  in  practice,  mutual  coupling  exists  between 
the  array  sensors.  Because  the  mutuals  change  the  sensor  impedances,  the 
gain  and  radiation  pattern  of  the  array  can  be  greatly  distorted.  In  sub¬ 
space  methods  this  can  significantly  alter  the  eigensystems  underlying  the 
estimation  procedures.  Gupta  and  Ksienski  [70]  investigated  the  effects  of 
mutual  coupling  on  the  performance  of  an  adaptive  array,  xu  their  treat¬ 
ment,  the  matrix  Zq  characterizing  the  mutual  coupling  between  the  sensors 
was  determined  using  a  mathematical  model  which  models  the  antenna  array 
consisting  of  m  sensors  as  an  (m+1)  terminal  network.  This  matrix  was  then 
used  to  determine  the  sensor  outputs  that  would  have  existed  had  there  been 
no  mutual  coupling.  A  compensation  scheme  was  then  developed  to  study  the 
method  known  as  beamforming.  Yeh  and  Leou  [71]  used  the  same  mathematical 
model  and  applied  it  to  the  MUSIC  algorithm.  Recall,  in  the  MUSIC  algo¬ 
rithm,  that  one  plots  the  inverse  of  the  correlation  between  the  noise  sub¬ 
space  En  and  the  directional  vector  a(0),  which  we  denote  by  ((En.a(0))2)- 
If  mutual  coupling  is  present,  Yeh  and  Leou  show  that  one  has  to  plot 
the  function  ((En.ZQ-^.a(0))2)~^.  Failure  to  do  so  results  in  severe 
degradation  of  the  estimates.  Shau  [44]  considered  the  case  of 
deterministic  signals  in  a  noise  free  environment  and  eliminated  the  ef¬ 
fects  of  mutual  coupling  for  the  method  known  as  MFBLP  (Modified  Forward 
Backward  Linear  Prediction).  In  this  chapter  we  develop  algorithms  to  ef¬ 
fectively  compensate  for  the  effects  of  mutual  coupling  when  using  the 
Matrix  Pencil  approach  with  ESPRIT  and  the  Moving  Window. 
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3.1  MODEL 

Consider  a  linear  array  of  m  dipoles  uniformly  spaced  at  a  dis¬ 
tance  D.  Each  dipole  is  of  length  t  and  has  a  radius  r  satisfying  the  con¬ 
dition  r«t.  A  load  is  attached  to  the  center  gap  of  each  dipole 
(Fig.  3-1).  Assume  there  are  d  narrowband  signals  impinging  on  the  array  as 
planar  wavefronts.  The  voltages  induced  by  the  assumed  signals  on  the  loads 
are  the  outputs  of  the  dipoles.  Induced  currents  will  appear  on  the 
dipoles.  These  currents  reradiate  and  generate  scattered  fields.  The  scat¬ 
tered  fields  then  induce  currents  on  the  neighboring  dipoles.  The  process 
of  induction  and  reradiation  causes  the  mutual  coupling  among  the  dipoles. 

Using  one  sinusoidal  expansion  and  weighting  function  per  dipole, 
the  method  of  moment  i  [45,46]  was  used  to  obtain  the  matrix  of  mutuals 
(Fig.  3-2).  Denote  the  current  distribution  by  J(z)  (assuming  longitudinal 
distribution  and  neglecting  all  other  distributions)  and  the  j-th  expansion 
function  by  fj(z).  Then 

m 

J(z)-E  I(j)f1(z)  (3.1-1) 

J-l 

where  I(j)  ;  j»l,  2,  .  .  . ,  m,  denotes  the  unknown  current  amplitude  to  be 
determined  on  each  dipole.  At  a  point  (y,z)  in  the  Y-Z  plane,  the  scattered 
field  is  given  by 

m 

Il<s>(y,z)-I  I(j)E(f>(y,z)  (3.1-2) 

j-l 
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where  Ej^s^(y,z)  is  the  scattered  field  from  the  j-th  dipole.  The  total 
field  will  then  be 

E(y,z)«E<inc)(y,z)  +  E<s>(y,z)  (3.1-3) 

where  E^nc^  is  the  incident  field.  Let  Ez  be  the  z-component  of  the  total 
field.  A  generalized  voltage  V(i)  induced  on  the  subsection  spanned  by  the 
function  fj(z)  can  be  defined  with  respect  to  a  weighting  function  wj(z)  as 
V(i)-F(Ez(y,z),wi(z))  (3.1-4) 

where  F  is  bilinear  with  respect  to  Ez(y,z)  and  wj(y,z).  Similarly,  we 
define  the  voltage  produced  by  the  incident  field  on  the  i-th  dipole  by 
V<inc)(i)-F(E^inc)(y,z),wi(z))  (3.1-5) 

and  the  voltage  produced  by  the  scattered  field  on  the  i-th  dipole  by 

V<s)(i)-F(E^s)(y,z),wi(z)).  (3.1-6) 

Thus,  the  total  voltage  introduced  in  the  i-th  dipole  is 
V(i)-V<inc>(i)+V<s)(i), 
which,  for  metallic  scatterers,  becomes 
V(i)-v(inc)(i)+V<s)(i)-0 
or 

V<inc>(i)«- V<s>(i).  (3.1-7) 

However , 


V<s>(i)  -F(  ^^I(j)E(j)(y,z),wi(z)) 


-  Z  I(j)F(E<f)(y,z),wi(z)). 

j-1 


The  total  impedance  between  the  i-th  and  j-th  dipoles  is  defined  to  be 
zi^-F(E(5)(y,z)fVi(z)).  (3.1-8) 

Thus, 
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•  f  Ill  • 


(3.1-9) 


V<s>(i)«E  -z^  I(j)  ;  i-1,2, .  . 

j-1 


In  matrix  notation 

V<s>  —  Z  I  (3.1-10) 

where 

v(s)T,[v<s)(1),v(s)(2),.  •  .,V^s>(m)] 

and 


It«[I(1),I(2),.  .  . ,I(m) ] . 

The  total  impedance  matrix  Z  can  be  decomposed  into  two  parts  as 

Z-Z0+ZL 

where 

Zg  is  the  generalized  impedance  matrix 

and 


Z^  is  the  load  matrix. 

Assuming  that  all  loads  are  loaded  with  the  same  load  zj.,  the  matrix  ZL  is 
given  by 


The  ij-th  element  of  Z,  therefore  is 

ziJ-Zij+ziSij , 

where  zjj  is  the  mutual  impedance  between  the  i-th  and  j-th  dipoles.  The 
total  voltages  induced  on  a  load  are  given  by 
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and 


V<L>  =  ZL  I 
I  -  Zl"1  v(L) . 

However, 

V(inc)  ,  zi  =  ZqZl_1V^l)  +  V<L) 

which  implies  that 

v<L>  -[i+ZqZl-1]-1  v<inc>.  (3.1-11) 

Let  H  be  the  matrix 

H«  [I+ZqZl"1].  (3.1-12) 

H  can  be  written  as 


H 


l+(z11/z1) 

(Z21/Zl> 


(z 12/Z1> 
l+(z22/zl) 


<zlm/zl> 

<z2m/zl> 


L 


<zml/zl) 


<zm2/zl> 


•  •  •  l+(zmn/zl)  J 


Thus,  when  incident  signals  are  impinging  on  the  array  and  in  the  presence 
of  additive  noise,  the  output  of  the  linear  array  will  be 

y<L>  *H_1  V<inc)  +  N. 

For  simplicity,  let 

X-V<inc> 


and 

Y«y<L). 

Ve  now  have  a  relationship  between  the  incident  signals  and  the  outputs  of 
the  array  which  is 

Y  -  H"1  X  +  N.  (3.1-13) 
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3.2  MINIMUM  MEAN-SQUARED  ERROR  ESTIMATION 

If  we  try  to  use  the  vector  Y  observed  at  the  output  of  the  array 
in  the  original  formulation  of  the  matrix  pencil  approach,  it  is  not  pos¬ 
sible  to  obtain  the  decomposition  needed  for  che  matrix  pencil.  An  estimate 
X  of  X  is,  therefore,  generated.  This  estimate  is  also  used  in  ESPRIT. 

Assuming  that  the  signals  and  noise  are  statistically  independent 
and  that  the  noise  components  are  uncorrelated  zero-mean  random  variables 
with  variance  a2,  the  minimum  mean-squared  error  linear  estimator  results 
when  the  error  (X-X)  is  orthogonal  to  the  observed  data  Y.  Let  X*  R  Y, 
where  R  is  to  be  determined.  Thus,  we  have 

E{(X-X).Yh}-0. 

However  X-  R  Y  which  implies  that 

E(X  YH}-RE{Y  f)  . 

Recall  that  Y  -  H*1  X  +  N.  Thus, 

E{X  YH}-E{X  (H-1  X  +  N)H}-E{X  ^(H*1)0 

and 

E{Y  Yh}»E{(H-1  X  +  NKH"1  X  +  N)0} 

*(H~^)E{X  X0)(H-1)0  ♦  a2  Im. 

Let  C  denote  the  correlation  matrix  of  Y.  Then 

C-E{Y  Y0}-  (H“1)E{X  X0}^-1)0  +  a2  Im. 

Thus, 

H  (C-o2  Im)-E{X  X0}^'1)0  -E{X  Y0}. 

Therefore, 

R-  E{X  Yh}(E{Y  Y0})'1 

R  -  H  (C-o2  Im)  C_1.  (3.2-1) 
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3.3  APPLICATION  TO  THE  MOVING  WINDOW 

With  the  new  formulation  of  the  moving  window,  pre-processing  the 
signals  received  at  the  output  of  the  array  will  allow  us  to  generate  the 
signals  that  would  have  resulted  had  there  been  no  mutual  coupling.  Recall 
that  in  the  presence  of  mutual  coupling,  the  received  signal  at  the  output 
of  the  array  can  be  modeled  as 

Y  «  H"1  X  +  N.  (3.3-1) 

X  represents  the  vector  of  incident  signals. 

Let  HI»H“1.  Consider  the  vector 

Y*  =  HI*  X*  +  N*  (3.3-2) 

where  *  denotes  complex  conjugate.  (L+l)  vectors  Xn  of  length  (m-L)  are 
then  formed  where 

In'*'  *  l  Yn  Y(n+1)  Y(n+2)  *  *  *  Y(n+m-L-l)  1* 

For  the  sake  of  clarity,  let  m-5  and  L*2.  Then 


yi*  hin*  hi12*  hi13*  hi14*  hi15*  1  r  XI*  1  r  nx* 

Y2*  hi21*  hl22*  hi23*  hl24*  hi25*  x2*  n2* 

y3*  ■  hisi*  hi32*  hi33*  hi 3^*  h^-35*  x3*  +  n3* 

y^*  ^^41*  ^*42*  hi43*  hi^^*  hi45*  X4*  04* 

Y5*  ■  ■  hi51*  hi52*  hi35*  hi54*  hi55*  -  •  x5*  -  •  n5* 

Note  that  this  matrix  equation  can  be  reformulated  as 

y3*’  hij3*  hii2*  hi33*‘  xj*"  D  0  hi^4*  x2*  0  0  hij3*  x3* 

y2*  *  hi2i*  hi22*  hi 23*  X2*  +00  hi24*  x3*  +00  hi25*  X4*  + 
y3*.  M31*  hi32*  ^^33*-  *3*-  &  ®  hi34*.  ■x4*-  -0  ®  hi33*.  ^5*- 

"hi2i*  0  O’  xj*"  hi22*  hi  23*  hi  24*  x2*  ®  0  hi  25*  x3* 

hi3j*  0  0  X2*  +  hi32*  hi33*  hi34*  x3*  +00  hi33*  X4*  + 

Jii41*  0  0.  x3*.  M42*  h*43*  ^44*.  *4*.  .0  0  hi43*.  .X5*. 
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V3*' 

hi31*  0  O’ 

*i* 

hi32*  0  O' 

r  *1 
x2 

M33* 

hi34* 

“35*' 

x3* 

r  *1 

n3 

?4* 

m 

hi41*  0  0 

x2* 

+ 

hi42*  0  0 

★ 

x3 

+ 

hi43* 

hi44* 

hi45* 

x4* 

+ 

★ 

n4 

ys*- 

M51*  0  0. 

*3*- 

hi52*  0  0. 

* 

Lx4  j 

^i53* 

hi54* 

hi55* 

■x5*- 

* 

Ln5 

Using  matrix  notation,  ve  have 

Yj  »  +  3^12  —2  +  Z3  +  Nj_, 

Y2  *  8^21  —1  +  ^22  —2  +  HI23  I3  +  —2’ 

Y3  *  HI31  +  HI32  Z2  +  HI23  Z3  +  N3. 

For  the  general  case,  it  can  be  shovn  that  Yjj  can  be  written  as 

(L+l) 

In  «  z  HIni  Zi  +  Nn  (3.3-3) 

i-1 


where 

IiT  *  I  xi  x(i+l)  x(i+2)  •  •  •  x(i+m-L-l)  1  *  2,  .  .  .,  (L+l) 

NiT  -  l  nt*  n(i+1)*  n(i+2)*  .  .  .  ]  ;  i-1,  2,  .  .  (L+l) 

and 


«ni 


“ni 


«ni 


hi 


ni 


hl(n+l)i 


ki(n+m-L-l)i 


hi 


ni 


hl(n+l)i 


0 

0 


hi 


(n+m-L-l)i  0 


0  0 
0  0 


•  Mn(i+ra-L-l)  ^ 

•  hi(n+l)(i+m-L-l) 


•  hl(n+m-L-l)(i+m-L-l)' 


.  0 

.  0 


.  0 


for  n>i 


0  0 


*  °  JjJn(i+m-L-l)  * 
0  hi(n+l)(i+m-L-l) 


hi 


(  n+m-L- 1 ) ( i +m-L- 1 ) 


for  n<i 


for  n=i 
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Note  that  HIn^  has  dimensions  (m-L)x(m-L).  Let  V  be  the  vector 

»  -  I  XlT  l2T  •  •  •  X(U1)T  1T- 

V  can  be  expressed  as 


[Mil 

“21 

hi12 

HI22 

*  *  *  HI1(L+1) 

•  *  •  **2(1+1) 

\h  1 

ll 

r  »i 
*2 

v  . 

• 

-  “(Ul)l 

•  • 

•  • 

®I(L+1)2  • 

•  •  * 

•  •  • 

•  •  '  HI(L+1)(L+1)  ■ 

■  i(ui)  ■ 

+ 

■  n(L+i)  ■ 

This  can  also  be  written  as 

V  «  %  Z  +  N'.  (3.3-4) 

Assuming  the  signals  and  noise  to  be  statistically  independent,  ve  get 
E[W  WH]  -  Hx  E[Z  Z11]^0  +  EJN'N'8].  (3.3-5) 

Therefore, 

E[Z  ZH]  -  (Hi)*1  (E[V  WH]-E[N'N'H])  ((H!)'1)11.  (3.3-6) 

However  E[Z  ZH]  can  be  expressed  as 


E[Z  ZH)» 


E(222iHI 

• 

E(2222h] 

•  • 

•  eI?iel!|] 

•  • 

ElEi5(L.i)!!l 

eIe2?(L.1)H) 

• 

EtglZl"]  „ 

•  • 

e(5i22HI  „  • 

•  •  • 

•  eIe(ui)Zlh)  eI?(l.1)5<ui>h1  j 

Thus,  after  partitioning  the  matrix  E[Z  ZHJ  into  a  total  of  (L+l)x(L+l) 
matrices,  each  having  order  (m-L)x((m-L),  the  matrices 


and 


L 

M  *  I  E[ Z±  ZiH  ] 
i-1 
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L 

N  -  r  B[Z(i+1)  ZiH  1. 

needed  in  the  formulation  of  the  matrix  pencil  (see  section  2.3.3)  are 
readily  obtained.  By  using  such  a  formulation,  we  can  effectively  com¬ 
pensate  for  the  effects  of  mutual  coupling  using  the  moving  window.  Recall 
that  the  i-th  incident  signal  is  given  by 

d 

Xi(t,0)-Z  sk(t)ai(9k);  i*l»  2,.  .  .,m 
k-1 


where 

ai(®k)  *  ak  eJ(i-1)*k 

♦k  *  -«Dsin(0k)/c  ;  k*l,  2,  .  .  .,d, 

to  is  the  center  frequency  of  the  plane  waves 

c  is  the  propagation  speed  of  the  waves 

D  is  the  sensor  spacing 

ak»a(0k)  is  the  beam  pattern  in  the  direction  of  the  k-th  emitter 

3.4  APPLICATION  TO  ESPRIT 

Three  different  arrangements  of  the  doublets  are  considered  in 
this  section. 

3.4.1  General  Array 

Ve  assume  that  we  have  2m  sensors  so  as  to  form  m  doublets.  Fur¬ 
ther,  we  assume  that  each  doublet  is  isolated  from  the  others  so  that 
mutual  coupling  exists  only  between  the  two  sensors  within  each  doublet 
(Fig.  3.3).  The  observed  signal  at  the  i-th  doublet  can  then  be  modeled  as 


vi 

wi 


xi 

■  yi  • 


'  nlj 
.  n2j 


i-1,  2, 


m. 


(3.4. 1-1) 
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Cone  mi  Case 
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Let  HIi-Hj-1.  HIj  can  be  written  as 

’  hilli  hil2i  ' 

.  hi21i  hi22i  J  . 


Then 

vj  «  hill^xj  +  hil2^yi  +nlj 

and  (3.4. 1-2) 

-  hi21jxj  +  hi22^yj  +n2^. 

Collecting  all  the  v^'s  in  a  vector  V  and  all  the  vj's  in  a  vector  V,  we 
have 

V  -  HIn  X  +  HI12  Y  +  Nx 

(3.4. 1-3) 

V  .  HI2i  X  +  hi22  Y  +  N2  , 


Where  HI^,  HI^2,  HI2j^,  HI22,  X,  Y,  and  N2  are  given  by 


HIn  -  diag  (hill!,  hill2,  .  .  hillm 

HI12  *  diag  (hil2! ,  hil22,  .  .  .,  hil2m  }, 

HI21  »  diag  (hi21lf  hi212,  .  .  .,  hi21m  }, 

HI22  -  diag  (hi221 ,  hi222>  .  .  hi22m  }, 

2  “  lxl»  x2»  *  •  •»  xml^» 

Y  -  [yi ,  y2,  .  .  xjT, 

5l  -  (nl^,  nl^,  .  .  .,  nlm]T, 

N2  *  { n2  ^  t  n2  ^ ,  •  •  • »  n2nj]'^. 

Consider  the  vector  Z  defined  by 

Z  -  [VT  VT  }T. 

Z  can  be  written  as 


'  Hill  Kl2  ' 

'  X  ' 

— 1  ' 

.  hi2i  hi22  . 

.  Y  . 

n2  . 

(3.4. 1-4) 
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Assuming  that  ‘he  signals  and  noise  are  statistically  independent  and  that 
the  noise  components  are  uncorrelated  from  sensor  to  sensor  with  covariance 
matrix  where  l2m  is  the  (2mx2m)  identity  matrix,  then  C22  =  E[Z  ZH, 

is  given  by 


' zz 


Hln  HIi2 
hi21  hi22 


E[X  XH]  E[X  YH] 
El Y  XH]  E[Y  YH] 


Hill  HI12 

H!2i  HI22  J 


H 


+  ff^I 


2m- 


(3.4. 1-5) 


Let  H2  be  the  matrix 


h2 


‘  Hill  HI12  ' 

■  HI21  HI22  J  • 


Then 


(H2)-!  <C22  -o2I2m)  ((H2)_1)H 


'  E[X  XH]  E[X  YH]  • 

.  E[ Y  XH]  E[Y  YH1  . 


(3.4.1 -6) 


Having  recovered  the  matrix  on  the  right  side  of  equation  (3.4. 1-6),  the 
matrices  M»E[X  XH]  and  N«E[X  YH]  can  be  identified.  Recall  that  incident 
signals  are  expressed  as 
d 

xj(t)  -  E  sk(t)gi(ejc) 
k-1 


d 

yt(t)  -  E  sk(t)e-J<wA/c>sin<0k> 
k-1 


where  g^S^)-  is  the  gain  response  of  the  i-;.h  --ensor  to  a  source  arriving 
at  angle  8^.  The  matrices  H  and  N  can  be  decomposed  as 

M  -  GSGH  and  N  -  GS*HgH  ,  (3.4. 1-7) 


5/ 


where  G,  S  and  f  are  given  by 


G  -  [gi  g2  •  •  •  Sd  1  *  Eain  matrix 
S-E[S  SH], 

ST  ■  [s^  S2  •  .  .  s<jl  -  impinging  signal  vector, 

♦  -  diag  [  eJ*l,  eJ*2,  .  .  e^d  ], 

4^— (ctA/OsinCd^)  ,  k=l,  2,  .  .  d 

Therefore,  the  effects  of  mutual  coupling  have  been  eliminated  and  the  rank, 
reducing  values  of  the  matrix  pencil  (M-XN)  are  given  by 

^  -  e-j<“A/c)',in<ek)  ;  k-1,2,  .  .  .,d.  (3. 4. 1-8) 

and  the  angles  of  arrival  of  the  sources  are  given  by 

9k  -  sin-1{jln(Xk)/(wA/c))  ;  k-1,2,  .  .  .,d.  (3.4. 1-9) 


3.4.2  Linear  Array 
3.4.2. 1  OVERLAPPING  CASE 

Consider  a  linear  array  of  (m+1)  seniors  and  assume  there  are  d 
(d<m)  narrowband  sources  located  at  angles  6^;  k-1,.  .  .,d.  In  this  case  we 
consider  two  sub-arrays  consisting  of  the  first  m  sensors  and  the  last  m 
s  isors  (Fig.  3.4a).  The  observed  signal  vector  at  the  output  of  the  array 
can  be  written  as 

Y  -  H*1  X  +  N.  (3.4.2. 1-1) 

Let  HI  can  be  written  as 


Thus ,  i f 


and 


‘  hin  hi12 

•  •  •  ^l(m+l) 

ai21  hl22 

•  •  *  hi2(m+l) 

■  hl(m+l)l  hi(m+l)2  * 

•  •  •  hi(m+l)(m+l)  • 

Ii-(yi  y2  •  •  •  ymlT. 
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ve  can  vrite 


I2-IY2  Y3  •  •  *  y(ra+l))T* 


11  m  HI11  *1  +  HI12  X2  +  Nx  (3.4.2. 1-2) 

and 

12  m  ’{I21  *1  +  aI22  X2  +  N2  ,  (3.4.2. 1-3) 

where  HIn»  Hl^t  HI2i»  HI22,  Nj^,  N2,  X^.and  X2  are  given  by 

Xi-[xi  x2  .  .  .  xm]T, 

X2*(x2  X3  .  .  .  X(m+l)l^* 

Hln  -  [hilli  hill?  .  .  .  hillm  ], 
hilli*  [hin  hi2i  .  .  .  himi)T;  i«l,...,m, 

hi12t  -[00...  0  hill(m+1)  ], 

HI21t  »  [hi221  0  .  .  .  0  ], 

HI22t  «  [hi22,  hi22?  .  .  .  hi22(m+1)  ], 

hi22j*  [ h±2i  h^3i  •  •  •  hi(m+l)ilT*  •  •  • » (m+1), 

5l  *  [nx,  n2,  .  .  .,  nm]T, 

T 

N2  -  [n2,  n3,  .  .  .,  n(m+1)]  . 

Consider  the  vector  Z  defined  as 

Z  -  HlT  I2T  1T- 

Z  can  be  written  as 


'  Hln  HI12  ' 

'Xi  ' 

*5i  ' 

■  HI21  hi22  - 

.  x2  . 

•52  ■ 

(3. 4. 2. 1-4) 


Assuming  that  the  signals  and  noise  are  statistically  independent  and  that 
the  noise  components  are  uncorrelated  from  sensor  to  sensor  with  variance 
<tZ,  then  Czz  *  E[Z  ZH]  is  given  by 


'zz 


rain  “12 
U21  hi22. 


nix^H]  eix^H] 
£[x2Xih1  e[x2x2h]J 


BI11  hi12 
jil2i  hi22 


-|H 


+  <7^ 


^m  ^m 


ImJ 


(3.4.2. 1-5) 
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Let  H3  and  [I]  be  the  matrices 


and 


■3 


’  HIU  HI12  ■ 
.  HI2i  HI22  ■ 


[I] 


*m  ^m 

-  -^m  ■  • 


where  Im  is  the  identity  matrix  and  Ilm  and  I2m  are 


0  0  0 
10  0 
0  10 


0  0  ‘ 
0  0 
0  0 


•  •  • 

0  0  0 

.000 


0  0 
1  0  . 


12. 


HmT 


Then 


(H3)-1  (Czz  -<*i])«i3)-l)H. 


'  ElXi  ?1H]  EtXi  X2H]  ' 

•  El?2  2lHl  Ef?2  ^2HJ  ■ 


(3.4.2. 1-6) 


Having  recovered  the  matrix  on  the  right  side  of  equation  (3. 4. 2. 1-6),  the 
matrices  M«E(XjXiH]  and  N»E[X^X2H]  can  be  identified.  Recall  that  the  i-th 
incident  signal  is  given  by 


d 

xj(t,9)»E  sjt(t)ai(9j,);  i-lt  2,.  .  .,(m+l) 
k-1 


where 


ai(®k)  *  ak  e^^"^^^k 

♦k  .  -wDsintOj^/c  ;  k-1,  2,  .  .  .,d, 
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(o  is  the  center  frequency  of  the  plane  waves 
c  is  the  propagation  speed  of  the  waves 
0  is  the  sensor  spacing 

ak“a(®k)  *s  beam  pattern  in  the  direction  of  the  k-th  emit¬ 
ter. 

It  can  be  shown  that  M  and  N  have  the  decompositions 

M  =  ASAH  and  N  =  AS^A3  (3.4.2. 1-7) 

where  A,  S  and  *  are  the  following  matrices 
S-E[S  SH], 

ST={si,  .  .  .  .sj}  impinging  signal  vector, 

A  “  lal  a2  *  *  1  id  1 

aA  -  [  a(9i)  a^eJ+i  .  .  .  a^e^i] , 

*  a  diag  (  eJ^l,  .  .  .,  e^d  ]. 

Therefore,  the  effects  of  mutual  coupling  have  been  eliminated  and  the  rank 
reducing  values  of  the  matrix  pencil  (M-AN)  are  given  by 

\  -  e-j(«A/c)sin(9k).  k,if2,  .  .  .,d.  (3.4. 2. 1-8) 

The  angles  of  arrival  of  the  sources  are  given  by 

8^  a  sin"1{jln(Ak)/(uA/c));  i-1,2,  .  .  .,d.  (3.4.2. 1-9) 

3. 4. 2. 2  NON-OVERLAPPING  CASE 

In  tuis  case  a  linear  array  composed  of  2m  sensors  is  used.  Two 
neighboring  sensors  are  considered  as  a  doublet  so  thtt  m  doublets  are 
formed  (Fig.  3.4b).  Let  there  be  d  (d<m)  sources.  Again,  the  received  sig¬ 
nal  at  the  output  of  the  array  is  modeled  as 

Y  «  HI  X  +  N  (3. 4. 2. 2-1) 

where  HI  is  given  by 


62 


■  hill 

hi12  .  . 

•  '  *)*l(2m) 

hi2i 

hi22  •  • 

•  *  hl2(2m) 

•  •  • 

■  hi(2m)l 

hi(2m)2  •  * 

•  •  • 

*  •  hi(2m)(2m)  ■ 

Let  vj  and  vx  be  the  signals  received  at  the  i-th  doublet. Then 


vi  *  y(2i-l) 

and  {  i*l *  2,  .  .  «  ,  m. 

wi  *  y<2i) 

Collecting  all  the  vj's  in  a  vector  V  and  all  the  v^'s  in  a  vector  V  ,  ve 
have 

V  »  Hi!!  Xx  +  HIi2  X2  +  Nx  (3.4. 2. 2-2) 

and 

W  -  HI2i  Xi  +  HI22  X2  +  N2  ,  (3. 4. 2. 2-3) 

where  Hi!!,  0I12»  HI21»  H*22*  ?1»  ?2»  and  I?2  are  8*ven  by 

2l-txl  x3  •  •  •  X(2m-1)]T» 

X2-[x2  x4  .  .  .  x(2m)]T, 

HInT  -  [ hi  1 1 1  hill?  .  .  .  hillm  ], 

hilli-  [hi(2i_i)i  hi(2i_1)3  .  .  .  hi(2i_i)(2m_i) ] ;  i«l,  2,  .  .  .,  m, 
hi12T  *  [hi!2i  hi!2?  .  .  .  hi!2m  ], 

hil2i-  [hi(2i_i)2  hi(2i_1)4  .  .  .  hi(2i_1)(2m) ] ;  i»l,  2,  .  .  .,  m, 
HI2iT  »  [ hi 2 1 t  hi212  .  .  .  hi21m  ], 

hi21im  Ihi(2i)l  hi(2i)3  •  •  •  hi(2i)(2m-l)l;  ial»  2,  .  .  .,  m, 

HI22T  -  [hi22i  hi22?  .  .  .  hi22m  J, 

hi22i»  [hi(2i)2  hi(2i)4  •  •  •  hi(2i)(2m)];  i«l,  2,  .  .  .,  m, 

Nx  -  (ni,  n3,  .  .  .,  n^2m_i j ]T, 

T 

N2  -  [n2,  n4,  .  .  . ,  n2m]  . 
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Consider  Che  vector  Z  defined  as 

Z  -  [VT  VT]T. 

Z  can  be  vritten  as 


■  HIU  HI12  • 

'll  ' 

’5l  ‘ 

z  - 

+ 

.  HI21  0I22  - 

■  x2  - 

■N2  - 

(3. 4. 2. 2-4) 


Assuming  that  the  signals  and  noise  are  statistically  independent  and  that 
the  noise  components  are  uncorrelated  from  sensor  to  sensor  with  covariance 
matrix  o^l2m  where  l2m  is  the  (2mx2m)  identity  matrix.  Then  C22  *  E[Z  ZH] 
is  given  by 


'zz 


HTn  HI12 
“21  hi22 


■E[XiXiH]  EIX^] 
“X2X iH]  E[X2X2H] 


Mil  “12 

“21  HI22  J 


H 


+  o^I 


2m 


(3. 4. 2. 2-5) 


Let  H4  be  the  matrix 


H4 


’  HI11  HI12  ‘ 

.  HI2i  HI22  *  * 


Then 


(H 4)-1  (C2Z  -o2I2m)  ((H4)"1)h 


'  E[X!  Xi0]  E[XX  X2H] 
.  E[X2  X!H]  E[X2  X2H] 


(3. 4. 2. 2-6) 


Having  recovered  the  matrix  on  the  right  side  of  equation  (3. 4. 2. 2-6),  the 
matrices  M«E[XjXi0]  and  N«E[XiX2H]  can  be  identified.  The  i-th  incident 
signal  can  be  vritten  as 


d 

Xj(t,9)-E  Sk< Oa^Ofc);  i-1,  2,.  .  .,2m 
“  k-1 
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where 

ai(®k.)  ■  ak 

^  .  -«Dsin(0jc)/c  ;  k«l,  2,  .  .  .,d, 

(a  is  the  center  frequency  of  the  plane  waves 
c  is  the  propagation  speed  of  the  waves 
D  is  the  sensor  spacing 

ajc»a(0jt)  is  the  beam  pattern  in  the  direction  of  the  k-th  emit¬ 
ter. 

M  and  N  have  the  following  decompositions 

M  -  ASAh  and  N  *  AS^A3  ,  (3. 4. 2. 2-7) 

where  A,  S  and  #  are  given  by 

A  ■  §2  *  •  •  id  1 

H  «‘(  aOi)  a(9i)ej2*i  .  .  .  3(6^  (2m-2)*i  ]  f 

S-E[S  SHJ, 

ST-{sj,  .  .  .  ,sd]  impinging  signal  vector, 

*  -  diag  [  eJ*l  e^2  .  .  .  e^d  ], 

Therefore,  the  effects  of  mutual  coupling  have  been  eliminated  and  the  rank 
reducing  values  of  the  matrix  pencil  (M-XN)  are  given  by 

\  -  e-j (wA/cJsinC©^) ,  k.1,2,  .  .  .,d.  (3. 4. 2. 2-8) 

The  angles  of  arrival  of  the  sources  are  thus 

6^  .  sin-1{jln(Xk)/(«6/c));  k-1, 2,  .  .  .,d.  (3. 4. 2. 2-9) 

3.5  Computer  Simulation 

The  scenario  used  for  this  simulation  consisted  of  two  incoherent 
sources  (d-2)  which  are  incident  on  a  linear  array  consisting  of  eight 
uniformly  spaced  half  wavelength  dipoles  (m-8).  The  sources  are  assumed  to 


65 


be  located  at  9]>16°  and  @2*24°  .  The  noise  was  simulated  to  be  white 
Gaussian  vith  zero-mean  and  unit  variance.  The  sensor  spacing  was  assumed 
to  be  half  wavelength  such  that  coD/c  -  n  .  The  load  impedance  was  taken  to 
be  the  complex  conjugate  of  the  self  impedance.  The  statistics  were  derived 
from  50  runs  where  100  snapshots  were  taken  in  each  run.  The  results  are 
shown  in  Fig.  3-5  to  3-16.  In  these  figures,  (1)  represents  the  moving 
window,  (2)  and  (3)  represent  ESPRIT  for  the  linear  case  when  overlapping 
and  non  overlapping  arrays  are  considered,  respectively,  and  (4)  cor¬ 
responds  to  ESPRIT  used  in  a  general  case.  Without  compensation,  note  that 
all  algorithms  fail  to  accurately  locate  the  two  sources  due  to  the  dis tor¬ 
sion  introduced  by  the  mutual  coupling.  With  the  compensating  schemes  de¬ 
veloped  here,  all  algorithms  identify  the  locations  of  the  two  sources  cor¬ 
rectly.  However,  ESPRIT  used  in  a  Linear  Overlapping  Case  performs  much 
better  than  the  remaining  algorithms.  This  is  due  to  a  larger  array  aper¬ 
ture.  However,  our  objective  was  not  to  perform  a  comparison  between  the 
different  algorithms  but  to  derive  effective  methods  to  compensate  for  the 
mutual  coupling  effects.  This  has  been  achieved  and  it  is  shown  that  com¬ 
pensation  of  the  mutuals  is  likely  to  be  a  necessity  if  acceptable  per¬ 
formance  is  to  be  obtained  in  practice.  In  the  different  figures  for  the 
mean-squared  error  and  the  variance,  the  y-axis  is  defined  as 
y-10  log10(.). 

Let  0^  be  an  estimate  of  9  obtained  at  the  k-th  run  (K  is  the  number  of 
runs).  The  sample  mean  (ME),  the  sample  variance  (Var)  and  the  mean-squared 
error  (MSE)  are  defined  respectively  as 

K 

ME( 0)  -  ( 1/K)  E  9^, 
k-1 
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K 

Var(0)  -  (1/K)  I  (0j,-ME(0))2, 
k*l 

K 

MSE(0)  -  (1/K)  I  <  0W-0)2 . 
k=l 
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SAMPLE  MEAN 


Without  Compensation:  ANGLE=16 


SNR  (dI3) 

( 1 ) -Moving 

window 

(2  J -ESPRIT 

:  Linear  Overlapping  C a s  e 

(3) -ESPRIT 

:  Linear  Non  Overlapping 

C  a  s  c 

(4) -ESPRIT 

:  General  Case. 

Fig.  3.5 

Sample  Mean  of  the  Angle 

Estimate  at  lb 

Without  Compensation  for 

the  Mutuals. 

(1) -Moving  Window 

(  2 )  -  ESPRIT :  Linear  Overlapping  Case 

(3)  -ESPRIT :  Linear  Non-Overlapping  Case 

(4 )  -  ESPRIT :  General  case. 


Fig.  5.6  Sample  Mean  of  the  Angle  Estimate  at  16° 

With  Compensation  for  the  Mutuals  When  Using 
a  Minimum  Mean-Squared  Error  Estimation. 


SAMPLE  MEAN 


Direct  Method:  ANGLE=16° 


SNR  (dB) 


( 1 )  -Mov ing  W indow 

(2)  -ESPRIT:  Linear  Overlapping  Case 

(3) -ESPRIT:  Linear  Non  Overlapping  Case 

(4)  -ESPRIT:  General  case. 


Fig.  3.7  Sample  Mean  of  the  Angle  Estimate  at  16° 

With  Compensation  for  the  Mutuals  When  Using 
the  Direct  Method. 
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Without  Compensation:  ANGLH=16 


( 1)  -Moving 

(2)  -ESPRIT 
(5) -ESPRIT 
(41 -ESPRIT 


W  i  iulow 

Linear  Overlapping  Case 
Linear  Non  Overlapping  Case 
General  Case. 


Fig.  3.8  Mean-Squared  lirror  of  the  Angle  estimate 

at  16°  Without  Compensation  for  the  Mutuals. 
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MSE  (dB) 


1 1 - * - » - * - * - 1 

5  10  15  20  25  30 


SNR  (dB) 


(1)  -Moving  Window 

(2)  -ESPRIT:  Linear  Overlapping  Case 

(3)  -ESPRIT:  Linear  Non  Overlapping  Case 

(4)  -ESPRIT:  General  Case. 


Fig.  3.9  Mean-Squared  Error  of  the  Angle  Estimate 
at  16°  With  Compensation  for  the  Mut’  i 1 s 
When  Using  a  Minimum  Mean-Squared  lii"  or 
E  s  t  i  ma  t i on . 
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MSE  (dB) 


Direct  Method:  ANGLE  =  16 


SNR  (dB) 


( 1 )  -Moving 

(2)  -ESPRIT 
(5) -ESPRIT 
(4)  -ESPRIT 


w indow 

Linear  Overlapping  Case 
Linear  Non  Overlapping  Case 
General  Case. 


Fig.  3.11)  Mean -Squared  Error  of  the  Angle  Estnuat 
at  16°  With  Compensation  for  the  Mutual 
When  Using  the  Direct  Method. 
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SAMPLE  MEAN 


Without  Compensation:  ANGLH=24 


( 1 )  -Moving  W  indow 

(2)  -ESPRIT:  Linear  Overlapping  Case 

(3)  -ESPRIT:  Linear  Non  Overlapping  Case 

(4)  -ESPRIT:  General  Case. 


Fig.  3.11  Sample  Mean  of  the  Angle  Estimate  at 
Without  Compensation  lor  the  Mutuals 


Minimum  Mean-Squared  Error  Estimation:  ANGLE=24° 


(1)  -Moving  Window 

(2)  -ESPRIT:  Linear  Overlapping  Case 

(3)  -ESPRIT:  Linear  Non  Ov  e  r 1  a pp  L ng  Ca s c 

(4)  -ESPRIT:  General  Case. 


Fig.  3.12  Sample  Mean  of  the  Angle  Estimate  at  24° 
With  Compensation  for  the  Mutuals  When 
Using  a  Minimum  Mean-Squared  Error 
Estimation. 
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SAMPLE  MEAN 


Without  Compensation:  ANGLE=24  ° 


:  < - - — - - 1 - 1 - - - ■  _ _i 

5  10  15  20  25  30 

SNR  (dB) 


(1)  -Moving  Window 

(2)  -ESPRIT:  Linear  Overlapping  Case 

(3) -ESPRIT:  Linear  Non  Overlapping  Case 

(4)  -ESPRIT:  General  Case. 


Fig.  3.14  Mean-Squared  Error  o  i  the  Angle  Estimate 
at  21°  Without  Compensation  lor  t  lie 
Mu  tua  1  s  . 
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MSE  (dB) 


Minimum  Mean-Squared  Error  Estimation:  ANGLE  -  24° 


( 1)  -Moving 

(2)  -ESPRIT 

(3)  -ESPRIT 

(4)  -ESPRIT 


IV  i  n  d  o  w 

Linear  Overlapping  Case 
Linear  Non  Overlapping  Case 
General  Case. 


Fig.  3.15  Mean - Squa red  Error  of  the  Angle  Estimate 
at  24°  With  Compensation  for  the  Mutuals 
When  Using  a  Minimum  Mean-Squared  Error 
Estimation. 
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Direct  Method:  ANGLE  =  24  ° 


(1)  -Moving  Window 

(2) -ESPRIT:  Linear  Overlapping  Case 

(3)  -ESPRIT:  Linear  Non  Overlapping  Case 

(4)  -ESPRIT:  General  ease. 


Fig.  3.16  Mean-Squared  Error  of  the  Angle  Estimate 
at  24°  With  Compensation  for  the  Mutuals 
When  Using  the  Direct  Method. 
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CHAPTER  4 

EXTENSIONS  TO  WIDEBAND  SIGNALS 

Different  methods  for  estimating  the  angles  of  arrival  of  wideband 
signals  can  be  developed  depending  upon  the  approach.  In  this  section  we 
devise  three  new  techniques  for  the  matrix  pencil. 

4.1  TRANSIENT  SIGNALS 

In  this  section  ve  model  each  source  as  a  sum  of  decaying  exponen¬ 
tials.  This  representation  is  appropriate  for  non  stationary  signals.  Con¬ 
sider  a  linear  array  which  consists  of  m  identical  wideband  sensors 
uniformly  spaced  at  a  distance  A.  Assume  there  are  d  broadband  sources  im¬ 
pinging  on  the  array  as  planar  wavefronts  and  emitting  signals  whose  com¬ 
plex  envelopes  are  denoted  by  sj^t).  The  signal  received  at  the  i-th  sensor 
can  be  expressed  as 

d 

xj(t)»E  3(0^)  s^t-T^)  +nj(t);  i»l,2,.  .  .,m  (4.1-1) 

k-1 

where  is  the  time  delay  that  source  k  takes  to  travel  from  the 
reference  point  to  the  i-th  sensor.  Taking  the  reference  as  the  first 
sensor,  can  be  written  as 

■Cik-a-lHA/cJsinOk)  (4.1-2) 

where  c  is  the  speed  of  propagation  of  the  waves.  Assume  the  k-th  source 
can  be  represented  by  a  sum  of  exponentials  having  natural  frequencies  p^; 
1-1,  2,  .  .  . ,  Thus,  Sk(t)  can  be  written  as 
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(4.1-3) 


M<k> 

sk(t)  -  £  b]_k  ePlk1  , 
1-1 


where  the  coefficients  are  assumed  to  be  random.  Therefore,  the 
received  signal  xj(t)  can  be  expressed  as 


d  M<k> 

xj(t)  -  Z  Z  Clk  e^i"1>^lk  +  n± ( t ) ;  i-1,2,.  .  .,m,  (4.1-4) 

k-1  1-1 


where 

♦lit  -  -( A/c)p^ksin(0jc) ,  (4.1-5) 

and 

cik  *  a(®k)^lkePlkt.  (4.1-6) 

Given  the  data  collected  at  the  output  of  the  array,  the  problem  is  to 
estimate  the  angles  of  arrival  of  the  sources.  From  the  above  data,  a 
matrix  pencil  is  generated.  It  is  shown  that  the  rank  reducing  values  of 
this  pencil  are  related  to  both  the  angles  of  arrival  and  the  natural  fre¬ 
quencies  of  the  sources.  The  natural  frequencies  of  the  sources  are  assumed 
to  be  unknown  at  the  receiver.  Therefore,  these  natural  frequencies  have 
first  to  be  estimated  and  then  be  used  to  solve  for  the  angles  of  arrival. 
Thus,  a  simultaneous  estimation  of  the  natural  frequencies  and  the  angles 
of  arrival  is  needed.  To  do  so,  the  first  sensor  is  followed  by  an  equally 
spaced  tapped  delay  line  consisting  of  m  taps  with  successive  delays  of  T 
seconds.  In  addition  the  received  signal  at  the  i-th  sensor  is  delayed  by 
an  amount  of(i-l)T;  i— 2 ,  3,  .  .  . ,  m,  (Fig.  4-1).  The  signal  at  the  output 
of  the  h-th  delay  following  the  first  sensor  is 

d  M<*> 

yh-xi(t-(h-1)T)-E  Z  clk  e<h-1)nk+n1(t-(h-l)T);  h-0, 1 , . . . ,  (m-1)  (4.1-7) 

k-1  l-l 
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where 


Ylk  “  "Plk^*  (4.1-8) 

The  signal  at  the  output  of  the  delay  connected  to  the  i-th  sensor  can  be 
expressed  as 

d  M<k> 

21«xi(t-(i-l)T).r  I  clke<i-1><*Lk+nk)+ni(t-(i-l)T);i=2,...,m.  (4.1-9) 
k-1  1*1 


Let 


*ij“(+ij+nj) 


and 

d  . 

M  -  E  M<k>  . 
k-1 


(4.1-10) 


(4.1-11) 


At  this  point  we  assume  that  all  sources  have  distinct  natural  frequencies. 
This  condition  is  relaxed  later  on.  With  the  knowledge  of  the  parameter  M, 
we  form  (m-L+1)  vectors  (m-L+1)  vectors  ^  and  (m-L+1)  vectors  Z,^  of 
length  L  where 


M  <  L  <  (ra-M) 


and 


Xn  *  Ixn» 

xn+l»  * 

•  *  *  Xn+L-1 

1T, 

n-1 , 

2,  .  . 

In  *  l^n' 

^n+1*  • 

*  * »  ^n+L-1 

1T, 

n«l , 

2,  .  . 

In  *  Izn» 

zn+l»  * 

•  •»  Zn+L-1 

1T- 

n-1, 

2,  .  . 

. ,  (m-L+1) 


It  can  be  shown  that  Xjj,  Yj^  and  Zj,  can  be  decomposed  into 
-  *l«(n-l)  C  ♦ 

Yjj  -  A2r<n_1)  C  +  NYjj, 


Zn  *  ASY^-1)  C  +  NZj,, 


(4.1-12) 

(4.1-13) 

(4.1-14) 


where  Al,  A2f  A3,  #,  T,  Y,  C,  NXn,  NYn ,  NZn  are  given  by 
Al»(  Allri  Al_2, 1  •  '  *  J' 
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A2-t  A2ifi  A22>1 

•  •  •  A2H(d)>(j  1, 

A3-1  A3lf !  A32,l 

.  .  .  A3M<d)fd  ]> 

Alifj  -  I  1  e*ij 

.  .  .eCL-D+ij  ]T, 

A2ifj  -  [  1  eYij 

.  .  .e<L-1)nj  ]T, 

A3itj  -  1  1  e*ij 

.  .  .eO-Dfij  ]T, 

♦»diag{  e^ll  .  . 

.  e*M<d>d}, 

T«diag{  eYll  .  . 

.  eYH<d)d}, 

Y-diag{  e*ll  .  . 

.  e*M<d)d}, 

C  *  [cji,  c2i,  • 

•  cM(d)d  JT, 

NXjj*  [  nxn ,  nxn+^ , . . 

..,nxn+L-l 

NYn-tnyn,nyn+1,.. 

••»nyn+L-l  J* 

NZn*[nzn,nzn+1, . . 

••*nzn+L-l  1* 

nxj,  nyj  and  nzj  are  introduced  here  for  simplicity  of  notation.  Actually, 
they  are  given  by 

nxj  -  ni(t), 

ny4  .  n1(t-(i-l)T), 

and 


nzi  -  n^t-U-DT). 


Note  that  nx^  ■  ny^  *  nz^.  Six  matrices  Mj_,  Nj,  P^,  Qj,  Rj  and  Sj  are 
formed  where 


’  t 

1 

t 

1 

'  t 

1 

t 

1 

Mi  - 

Si  • 

•  •  ?(m-L) 

;  Ni  - 

x2  . 

•  •  S(m-Ul) 

1 

.  ± 

1 

1 

1 

.  A 

1 

i 
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Mi  can  be  rewritten  as 


Hj  *  A1  [  C  K  .  .  .  *<n»-L-l)  c  ]  +  (  NXj  NX2  .  .  .  NX(m_L)  ]• 


Similarly,  it  can  be  shown  that 

Nx  -  Al*  l  C  *C  .  .  .  C  ]  +  [  NX2  NX3  .  .  .  NX(m_L+l) 

Px  -  A2  [  C  rc  .  .  .  r<m-L-1)  C  ]  +  [  NYX  NY2  .  •  .  NY(m_L)  ], 

Qi  -  A2r  [  C  rc  .  .  .  r<m-L-1>  C  ]  +  [  NY2  NY3  •  •  .  NY(m-Ul)  ]' 

Rx  .  A3  I  C  1C  .  .  .  c  I  +  [  NZt  NZ9  .  .  .  NZ(m_L)  1» 

Si  -  A3T  I  C  tc  .  .  .  y(®-l-1)  c  ]  +  f  NZ2  •  •  •  NZ(m-Ul)  1- 

These  matrices  have  the  following  decompositions 

Mx  -  Al  C  U1  +  Nl'  and  Nx  -  Al  C  *  U1  +  Nl"  ,  (4.1-15) 

Pi  -  A2  C  U2  +  N2'  and  Qi  -  A2  C  T  U2  +  N2"  ,  (4.1-16) 

Rl  -  A3  C  U3  +  N3'  and  Sx  -  A3  C  Y  U3  +  N3"  (4.1-17) 

where  Ul,  U2,  U3,  C,  Nl',  N2',  N3',  Nl",  N2"  and  N3"  are  given  by 

UlT-[  Ullrl  Ul2,l  •  ■  •  UlH<d),d  I. 
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U2t-(  U2lfl  U22fl  .  .  .  U2„<d),d  ], 

U3T-l  U314  U32>1  .  .  .  U3„(d)>d  I, 

Uli,j  -  l  1  e*ij  •  •  .e<m-L-1)^ij  ]T, 

U2ifj  -  l  1  eYij  .  .  .e(m-L-l)rij  f , 
y3i?j  «  [  1  e*ij  .  .  .e(m-L-l)*ij  ]Tf 
C  -  diag  [cn,  c2i,  .  .  cM(d)d  ], 

Nl'-lNXi  NX2  .  .  .  NX(m_L)  ], 

N2»-[OT1  NY2  •  •  •  NY(m-L)  1* 

N3'»[NZi  N22  .  .  .  NZ(m_L)  ], 

N1"-(NX2  NX3  .  .  .  NX(m_ui)  1, 

N2"-(NY2  NY3  .  .  .  NY(m_L+1)  }, 

N3MNZ2  NZ3  .  .  .  NZ(m.ul)  ]. 

Assuming  the  signals  and  noise  to  be  statistically  independent  and  that  the 
noise  components  are  uncorrelated  from  sensor  to  sensor,  we  get 
ElM^MiJ-Ul0  VI  U1  +  L  a2  I(m_L)  , 

ElN^MiJ-Ul0  #  VI  U1  +  L  a2  I1(m_L)  , 

E(P1HP11-U2h  V2  U2  +  L  a2  I(m_L)  , 

BIP1HQ11-U2H  I*  V2  U2  +  L  o2  I1(m_L)  , 

ElR^l-ltt0  V3  U3  ♦  L  a2  I(m_L)  , 

E[S1®Ri]*U30  T®  V3  U3  ♦  L  a2  I1(m_L)  , 
where  I(m_L)  is  the  (m-L)x(m-L)  identity  me.trix  and  Ii(m_L)»  VI »  V2  and  V3 
are  the  matrices 


0  10  0 
0  0  10 
0  0  0  1 


Il(m-L)  “ 


0  0  0  0 

.0  0  0  0 


.  0 
.  c 
.  0 


.  1 

.  0  . 
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V1-E(CHA1HA1C J , 

V2-B{CHA2HA2C], 

V3-E[CHA3HA3C1. 

Now,  let  M,  N,  P,  Q,  R  and  S  be  the  matrices 

M  -  E[M1HM1 1  -  La2  I(m_L)  -  U1H  VI  U1  ,  (4.1-18) 

N  -  E[N1HM1 1  -  La2  VI  U1  ,  (4.1-19) 

P  -  E( Pin?! 1  -  L a2  I(m_L)  -  U2H  V2  U2  ,  (4.1-20) 

Q  -  ElQjHPi]  -  La2  Ii(m_L)  -U2H  TH  V2  U2  ,  (4.1-21) 

R  -  E{R1HR1)  -  La2  I(m_L)  -  U3H  V3  U3  ,  (4.1-22) 

S  -  E[ S1hR1 1  -  La2  Ii(m_L)  «U3H  t8  V3  U3  .  (4.1-23) 

Consider,  the  following  three  pencil  matrices  (M-XN),  (P-tiQ)  and  (R-vS). 

Note  that 

(M-XN)-(U1H  VI  U1)-X(U1H  VI  U1)-U1H(I-X*H)V1  U1  ,  (4.1-24) 

(P-nQ)-(U2H  V2  U2)-n(U2H  rH  V2  U2)-U2H(I-qrH)V2  U2  ,  (4.1-25) 

(R-vS)-(U3H  V3  U3)-v(U3H  V3  U3)-U3H(I-v'fH)V3  U3  .  (4.1-26) 


The  matrices  Ul,  U2,  U3,  ♦,  T,  and  ¥  are  all  of  rank  M  as  long  as  all  the 
are  distinct  and  L>M.  Defining 

L 

Epq,rs  “  E  pq_^rs^  ’ 

i  « 1 

L 

Gpq,rs  *  1  «xp{(i-l)(r*pq-rrs)}, 

km.  iml 

L 

^pq , rs  ”  E  ®xP((i-l)(¥  pq_¥rg)) 

and 

cpq,rs  “  Elcpqcrs  )» 

the  matrices  VI,  V2  and  V3  can  be  written  as 
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cufll  Fll,  11 

c12,ll  p12,ll 


LC„(d)M(d)fll  F„(d)M(d)tll 


Cll,M<d>M<d>  PllfM(d)M<d> 
c12,M<d>M<d>  F12,M(d>M<d> 


CM<d)M(d>,M<d>M<d)  FM(d)M(d)fM(d)M(d)  J 


V2- 


cll,ll  Gll,ll 
c12,ll  g12, 11 


L  CM(d)M(d)fll  GM(d)M(d)fll 


cll,M<d>M(d>  Gll,M<d>M(d> 
c12,M<d>M^d)  G12fM(d>M<d> 


CM<d>M<d>,M(d>M(d>  GM(d)H(d)fK(^)H(d)  J 


C11 , 11  H11,H 
c12, 11  H12 ,11 


clltM(d)M<d>  Hll,M(d>M(d> 
c12,M<d>H<d>  H12,M<d>M<d> 


l  C„(d)M(d)fll  H„(d)M(d)fll 


cM<d>M<d>,M(d>M<d>  HM(d)M(d)fM(d)„(d)  J 


It  is  easy  to  see  that  the  matrices  VI,  V2  and  V3  are  of  rank  M  even  in  the 
presence  of  fully  correlated  sources.  The  rank  of  the  pencil  (M-AN)  is 
decreased  by  1  whenever 

Aij-exp{-*ij*}-exp{p1j*(A/c)sin(9j)} ,  (4.1-27) 

for  i-1,  2,  .  .  . ,  M<k>, 
j , k-1 ,  2 ,  .  .  ■ .  d . 

The  rank  of  the  pencil  (P-hQ)  is  decreased  by  1  whenever 

hij-exp{-Yij*}-exp{pij*T))  ,  (4.1-28) 

for  i-1,  2,  .  .  . ,  M<k\ 

J  f  ic*  1>  2  f  •  •  *9  cl* 

The  rank  of  the  pencil  (R-vS)  is  decreased  by  1  whenever 

\»1j-exp{-^1j*)-exp{pij*(  A/c)sin(9j)}exp{pij*T} ,  (4.1-29) 

for  f-1,  2 . M<k>, 

j , k-1 ,  2 ,  .  .  . ,  d . 

Note  that  the  first  set  of  generalized  eigenvalues  gives  us  a  set  of 
coupled  values  of  natural  frequencies  and  angles  of  arrival.  The  second  set 
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solves  for  the  natural  frequencies.  With  these  two  sets  one  might  think 
that  the  problem  is  solved.  However,  there  is  some  ambiguity  in  choosing 
which  natural  frequency  goes  with  which  angle  of  arrival.  The  third  set 
solves  this  ambiguity  since  we  can  see  that  the  set  of  these  generalized 
eigenvalues  is  the  product  of  the  first  and  the  second  ;i.e; 

vij  ■  •  *Uj*  (4.1-30) 

Therefore,  the  ambiguity  is  removed  by  constructing  a  table  of  all  the  pro¬ 
ducts  of  Xij  and  r^.  These  products  are  then  compared  with  the  values  vrs. 
Once  a  product  is  matched,  that  natural  frequency  is  used  to  determine  the 
angle  of  arrival  of  the  source.  In  practice,  due  to  numerical  round  off  er¬ 
rors,  a  range  of  uncertainty  remains  since  the  products  and  the  rank  reduc¬ 
ing  values  of  the  third  set  do  not  match  exactly. 

In  the  above  algorithm,  note  that  knowledge  of  the  number  of  natu¬ 
ral  frequencies  is  important.  If  sources  have  common  natural  frequencies, 
then  these  frequencies  would  be  counted  only  once.  This  means  that  the  to¬ 
tal  number  of  distinct  natural  frequencies  M  should  be  replaced  by  a 
suitably  reduced  number. 

4.1.1  COMPUTER  SIMULATION 

In  this  simulation,  a  linear  uniformly  spaced  array  consisting  of 
12  sensors  was  used.  Two  sources  were  assumed  to  be  present  and  located  at 
angles  9i-16°  and  02-24°.  Two  cases  were  studied. 

Case  1. 

Source  |  Angle  of  arrival  |  Natural  frequencies  I 


1  1 

16° 

|  0.25-jO. 90 

;  0. 25+ jO. 90 

2  1 

24° 

|  0. 12-j0. 79 

;  0. 12+j0. 79 
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Source 

Case  2. 

|  Angle  of  arrival 

|  Natural  frequemcies 

1 

|  16° 

|  0. 25-j0. 90  ;  0.25+j0.90 

2 

|  24° 

|  0.25-jO. 90  ;  0. 25+j0. 90 

The  sensor  spacing  was  selected  such  that  A/c*Jl/0.90.  The  coefficients  cj^ 
were  assumed  to  be  independent  random  valuables.  The  additive  noise  vas 
generated  as  white  Gaussian  with  zero  mean  and  unit  variance.  100  snapshots 
were  considered  in  each  of  the  50  runs  simulated.  The  results  of  the  angle 
estimates  are  shown  in  figures  4.3  to  4.8.  It  is  clear  from  these  figures 
that  the  estimates  obtained  from  the  second  case  (assuming  2  common  poles) 
are  much  better  than  case  (1).  This  is  mainly  due  to  the  fact  that  the  nat¬ 
ural  frequencies  estimates  are  less  biased  in  this  case.  The  length  of  the 
window  is  smaller  which  results  in  a  better  noise  reduction  through  the 
singular  value  decomposition  (SVD).  Tables  4.1  to  4.6  give  the  poles 
estimates  with  their  variances.  It  is  apparent  here  that  due  to  the  fact 
that  fewer  poles  had  to  be  estimated  in  the  first  case,  the  estimation  pro¬ 
cedure  achieved  a  better  performance  than  in  the  second  case. 
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SNR  (dU) 


(1) -No  Common  Poles 

(2) -2  Common  Poles 


Fig.  4.3.  Sample  Mean  of  Angle  at  10 
Transient  Signals. 


SAMPLE  MEAN 


ANGLE=24 


(1) -No  Common  Poles 

(2) -  2  Common  Poles 


Fig.  4.6.  Sample  Mean  of  Angle  at  24° 
Transient  Signals. 
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SNR  (dB) 

j  Sample  Mean 

Sample  Variance 

30 

|  0.2500616-jO. 9001127 

7 . 8816270e-4 

25 

|  0.2501127-j0. 9002026 

1 . 4029399e-3 

20 

|  0.2502106-j0. 9003684 

2 . 5036666e-3 

15 

|  0.2504140-j0. 9006952 

4. 5110551e-3 

10 

|  0.2509581-j0. 9014653 

1 .9616835e-2 

5 

|  0. 2519363-j0. 9021969 

1 . 9616835e-2 

0 

|  0.2310321-j0.969ul42 

0.1542990 

Table  4..1 

Sample  Mean  and  Variance  of  Che  pole  P| j =0 . 25- j 0 . 90 . 

(No  Common  poles) 


SNR  (dB) 

|  Sample  Mean 

1 

Sample  Variance 

30 

|  0. 2499895-jO. 8999775 

1 

8 . 9807872e-4 

25 

|  0.2499841-j0. 8999586 

1 

1 . 6006386e-3 

20 

j  0.2499823+j0. 8999237 

! 

2 . 6825261e-3 

15 

|  0. 2500210- jO. 8998682 

I 

5 . 1696543e-3 

10 

|  0.2503774-j0. 8998694 

1 

9 . 6555240e-3 

5 

|  0.2538730-j0. 9004380 

1 

2 . 2066321e-2 

0 

|  0. 2310484+j0. 9300936 

1 

9 . 0392388e-2 

Table  4.2 

Sample  Mean  and  Variance  of  the  pole  pi2=0- 25-jO. 90. 

(No  Common  poles) 
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0. 2179931- jO . 8058360 


•  Table  4.3 

Sample  Mean  and  Variance  of  the  pole  pn=Q.  12-j0. 79. 

(No  Common  poles) 


SNR  (dB) 

|  Sample  Mean 

Sample  Variance 

30 

|  0.1198654+j0. 7897864 

1.6900700e-3 

25 

|  0. 11977 49 +j 0.7896364 

3 . 0200828e-3 

20 

|  0.1196503+j0. 7894068 

5 . 4612877e-3 

15 

|  0.1195874+j0. 7891387 

1.0212054e-2 

10 

|  0. 1204553+j0. 7893061 

2 . 0683207e-2 

5 

|  0. 1317147 +j 0.7937539 

5.0145626e-2 

0 

|  0.2189744+j0. 8155535 

0.1081801 

Table  4.4 

Sample  Mean  and  Variance  of  the  pole  Pn=Q-  12-j0. 79 . 

(No  Common  poles) 
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SNR  (dB) 

|  Sample  Mean 

i 

Sample  Variance 

30 

|  0. 2500181 -jO. 8999961 

1 

25 

|  0.2500326-jG. 89Q993i 

I 

2 . 7576°82e-4 

20 

|  0. 2500592- jO . 8QQQ88 1 

1 

4.9048965e-4 

15 

|  0.2501090-jO. 899979? 

1 

8.  7285927e-4 

10 

|  0. 2502039- jO. 8999674 

1 

1 . 5552028e  -  3 

5 

|  0. 2504041-J0. 8999538 

I 

2  - 7789618e-3 

0 

|  0. 2508400- jO. 0999541 

1 

4.9983007c  3 

Table  4.3 

Sample  Mean  and  Variance  of  the  pole  pj^O.  25- }0.90. 

(2  Common  poles) 


SNR  (dB) 

|  Sample  Mean 

|  Sample  Variance 

30 

|  0. 25002 19* JO. 9000 132 

|  1 . 50806 2 9e- 4 

25 

|  0. 2500393* jO. 9000236 

|  2 . 68 1 307  7e-4 

20 

|  0. 2500712+ jO. 9000421 

|  4 . 7  6  5  4  8  7 1  e  -  4 

15 

|  0. 2501 304 +j 0.90007 52 

|  8 . 4  6  7  4  1 9  3  e  -  4 

10 

|  0. 2502445+ jO. 9002433 

|  1 .50389.\6e-3 

5 

|  0.2504744+jO. 9002433 

|  2 . 6 707326c- 3 

0 

|  0. 2509 71 3+ jO. 900444 7 

|  4 . 749 i ;r oP  i 

Table  4 .  <> 

''•ample  Mean  and  Variance  of  the  pole  pj  j -0. 23*  jO.'MK 

(2  Common  poles) 
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4.2  VTDB  SBNSB  STATIONARY  SIGNALS 

Consider  the  same  configuration  as  in  section  4.1.  Let  m  be  the 
number  of  videband  sensors  in  the  linear  uniformly  spaced  array  and  A  be 
the  sensor  spacing.  Assume  there  are  d  (d<m)  videband  sources  located  in 
the  far  field  so  that  planar  waves  arrive  at  the  array.  The  sources  s^(k) 
are  modeled  as  the  stationary  output  of  a  finite  dimensional  linear  system 
driven  by  white  noise  sequences  e^(k);  1-1,  2,  ...,  d,  and  k  is  a  discrete 
time  index  (53,54|.  Denote  the  transfer  function  of  the  1-th  linear  system 
as  h}(z).  Let  the  spatial  array  be  modeled  in  terms  of  the  impulse  response 
of  each  element  in  the  array.  The  array  response  to  a  unit  impulse  arriving 
at  the  array  from  direction  6  will  then  be  represented  as  the  impulse 
response  of  this  system.  The  combination  for  the  1-th  source  can  therefore 
be  modeled  as  a  single  system,  a^(z),  driven  by  a  vhite  noise  source  se¬ 
quence  (Fig.  4.2).  Let  gj(k,9^)  be  the  response  of  the  i-th  sensor  to  a 
unit  impulse  coming  from  direction  9^.  The  received  signal  at  the  i-th 
sensor  can  then  be  written  as 

d 

X}(k)-£  g1(k,91)*s1(k)  ♦  n t < k) ;  i-1,  2 . m  (4.2-1) 

1-1 

d 

xj(k)-I  gjik, 9})*hi(k)*ej(k)  +  nj(k);  i-1,  2,  .  .  .,  m  (4.2-2) 
1-1 

where  nj(k)  is  the  additive  noise  assumed  to  be  uncorrelated  with  the  emit¬ 
ter  signals  and  *  denote  the  operation  of  convolution.  Define  a^(k,9^)  as 

ai<1>(k,91)  -  gi(k,91)*h1(k).  (4.2.3) 

Let  A^)(z)  be  the  z-transform  of  a^^(k).  Assume  this  transfer 
function  to  be  a  rational  function  with  the  degree  of  the  numerator  being 
less  than  the  degree  of  the  denominator.  It  can  be  expressed  as 
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White  Noise 
Sequence 
eL  (k) 


H  L  (  2  ) 
hjtk) 


1-th  Source  Signal 
Sj^  (k)  »h  L  (k)  *c  L  (k) 


Fig.  4.2a  1-th  emitter  source 


i-th  sensor 
output 


Fig.  4.2b  Output  of  i-th  sensor  to  a  source 
coming  fron  direction  0^ 


White  noise 
sequence 


ex(k) 


AiU.V 

a^k.Gj) 


i-th  sensor  output 


Fig.  4.2c  Combination  of  the  two  systems. 


(4.2-4) 


M<1> 

A<1)(z)  -  I  h^/d-p^z-1) 
r-1 

where  p^r  is  a  complex  number  with  magnitude  less  than  1.  Note  that  these 
poles  do  not  belong  to  h^(k)  nor  to  g^(k,9^).  They  are  however,  a  mixture 
of  both.  Ve  will  loosely  refer  to  them  as  poles  of  the  sources.  Thus, 
a<L>(k)  can  be  written  as 

M<1> 

a<1>(k)  -  Z  hlr  <Pir)k.  (4.2-5) 

r-1 

However,  aj^^(k,9^)  is  related  to  a^^(k)  through 

a^1)^,^)  -  •<1)[k-(i-l)(A/cTs)sln(e1))  (4.2-6) 

where  Ts  is  the  sampling  period  and  c  is  the  propagation  velocity  of  the 
plane  waves.  Substituting  for  the  expression  of  a^^(k),  aj^^k,©^)  can 
then  be  written  as 

M<1> 

ai<1>(k,01)  -  Z  hlr  (plr)k  (Plr)"(i"1)(A/cTs)sin(el).  (4.2-7) 
r-1 


Let  be 

♦lr  -  (Pir)-(A/cTs>sin<0l>. 
aj(k,©2)  can  then  be  rewritten  as 


M*1) 

a^k,^)  -  Z  hlr  (plr)k  (♦lr)(i-1>. 
r-1 


(4.2-8) 


(4.2-9) 


Assuming  an  N  point  sequence,  the  received  signal  at  the  i-th  sensor  can  be 
expressed  as 
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(4.2-10) 


d  M<1) 

Xl(k)  -  E  E  hlr  (♦lr)<i~1>  <plr)k  *  e^(k)  +  ni(k) 
1-1  r-1 

1*1 )  2 ,  .  *  • ,  m 
k-0,  1,  .  .  ( N— 1 ) 


Let  Xj(n)  be  the  discrete  Fourier  transform  of  Xj(k).  This  is  given  as 


Xj(n)  *  E  xj(k)  (2n/N)nk> 
k-0 

N-l  (  d  M*1) 


(4.2-11) 


N-l  (  d  ) 

Xt(n)  -EE  E  hlr  (♦lr)<i'1>  (plr)k  *  ex(k)  e-J<2n/N)nk 
k-0  l  1-1  r-1  ) 


(4.2-12) 


+  E  nt(k)  e-j(2,l/N>nk. 
k-0 


i-1 ,  2  j  .  .  . ,  m 

n-0,  1,  .  .  . ,  (N-l) 


Let  DFT{ . )  denote  the  discrete  Fourier  transform  operator.  Then,  by  defini¬ 


tion, 


Hir(n)  -  hlr  E  <plr)k  e-j<2n/N)nk  «  hlr  DFT{(plr)k), 
k-1 


Sx(n)  -  DFT{e1(k)) 


Ni(n)  .  DFT{ni(k)) 


It  follovs  that 


d  M*1) 

Xt(n)  -  E  E  Hlr(n)  (♦lr)<i-1>  S^n)  ♦  N*(n) 
1-1  r-1 

i-1,  2,  . 
n-0 ,  1 ,  . 


(4.2-13) 


. ,  tn 

.,  (N-l), 


Given  this  set  of  data,  the  objective  is  to  solve  for  the  angles  of  arrival 
of  the  sources.  It  can  be  shovn  that  the  rank  reducing  values  of  the  matrix 
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pencil  generated  from  this  data  are  functions  of  both  the  angles  of  arrival 
and  the  poles  of  the  sources.  However,  the  poles  of  the  sources  are  assumed 
to  be  unknown  at  the  receiver.  Thus,  these  poles  have  first  to  be  estimated 
and  then  be  used  to  solve  for  the  angles  of  arrival.  That  is  the  reason  why 
the  same  configuration  as  in  section  4.1  can  be  used  to  solve  this  problem. 
Therefore,  the  first  sensor  is  followed  by  an  equally  spaced  tapped  delay 
line  consisting  of  m  taps  with  successive  delays  of  T  seconds.  In  addition 
the  received  signal  at  the  i-th  sensor  is  delayed  by  an  amount  of(i-l)T; 
i-2,  3,  .  .  . ,  m  (Fig.  4.1).  The  signal  at  the  output  of  the  h-th  delay 
following  the  first  sensor  is 

yh(k)-x1(k-(h-l)T) 
d  M<1) 

-  E  E  hlr  (Plr)k  (plr)*<h-1>T  *  e1(k-(h-l)T)  +  ni(k-(h-l)T)  (4.2-14) 
1-1  r-1 

d  M<1> 

-  E  E  hlr  (Plr)k  (nr)(h-1>  *  eX(k-(h-l)T)  ♦  ni(k-(h-l)T)  (4.2-15) 
1-1  r-1 

where 

TLk  *  "PlkT*  (4.2-16) 

The  signal  at  the  output  of  the  delay  connected  to  the  i-th  sensor  can  be 
expressed  as 

zs(k).xs(k-(s-l)T) 
d  M<k> 

-  E  E  hlr(plr)k  (^  Yir)<S_1>  *  e1(k-(s-l)T)  ♦  ns(k-(s-l)T)  (4.2-17) 
1-1  r-1 

Let 

*lr-<*lr><rir>*  (4.2-18) 

It  can  be  shown  that  the  discrete  Fourier  transforms  Y^(n)  and  Zs(n)  of 
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yjj(k)  and  zs(k)  are  respectively 


d  M<1> 

yh(n)  -  I  £  Hlr(n)  (Yir)^-1^  Sx'  (n)  +  Nh'(n)  (4.2-19) 

1-1  r-1 


Zs(n) 


h*  If  2  I  a  a  •  ,  nt 

n-0,  1,  .  .  (N-l ) 


d 

£ 

1-1 


M<1) 

£  Hlr(n)  (^lr  rir)^s_1)  S^n)  +  Nh"(n) 
r-1 

s— 1 ,  2,  .  . 
n-0,  1,  .  . 


(4.2-20) 


m 

(N-l) 


Assuming  that  the  sources  do  not  share  any  common  pole,  let  H  be 


d 

M  -  £  M*1)  .  (4.2-18) 

1-1 


As  in  section  4.1,  given  M,  ve  form  (m-L+1)  vectors  Xn,  (m  -L+l)  vectors  Yj, 
and  (m-L+1)  vectors  of  length  L  where 

M  <  L  <  (m-M) 


Xv(n)  -  [xv(n),  xv+1(n),  .  .  .,  xv+L_1(n)  ]T,  v-1,  2, 

Yv(n)  -  iyv<n>»  yv+i(n>»  •  •  •»  yv+L-i<n>  lT»  v-1*  2» 

Zv(n)  -  lzv(n),  zv+1(n),  .  .  .,  *v+L-l(n>  1T-  Vml>  2» 

It  can  be  shown  that  Xv,  Yv  and  Zv  can  be  decomposed  into 

X^n)  -  A1  H(n)  ♦(v-1)  S(n)  +  NXv(n), 

Yv(n)  -  A2  H(n)  T^"1)  S'(n)  +  NYv(n) , 

Zv(n)  -  A3  H(n)  Y<v-1>  S"(n)  +  NZv(n), 
where  Al,  A2,  A3,  ♦,  T,  Y,  C,  NXV,  NYV,  NZV  are  given  by 
Al-(  Allfl  A12>j  .  .  •  AlM^d>,d  J' 

A2*{  A2i,i  A22,l  *  •  •  ^M^.d  ], 

A3— l  A3iti  A32 , 1  *  •  •  Mh(d>,d 


(m-L+1) 

(m-L+1) 

(m-L+1). 


(4.1-19) 

(4.1-20) 

(4.1-21) 
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Mij  -  I  i  ♦ij  •  •  .  (♦ij)L-1  ]T, 

-  (  1  rij  •  •  •  (rij )L_1  ]T, 

A3ij  -  I  1  .  •  .<*ij)L-1  ]T, 

H(n)-diag{  Hlfl(n)  H1>2(n)  .  .  .  Hd>M(d)}, 

*-diag{  .  .  .  <frdM(d)}, 
r-diag{  "^11  .  .  .  Y(jM(d)} , 
f-diag{  .  *dH(d)}, 

S(n)  -  [S^n)  .  .  .  S^n).  .  .  Sd(n)  .  .  .  Sd(n)  ]T, 

S'(n)  -  tSj'Cn)  .  .  .  Sx'Cn).  .  .  Sd'(n)  .  .  .  Sd'(n)  ]T, 

S"(n)  -  [S1"(n)  .  .  .  Si* (n) .  .  .  Sd"(n)  .  .  .  Sd"(n)  JT, 

NXv(n)-[Nxv(n),  Nxv+1(n),  .  .  ,  Nxv+L_1(n)  ], 

NYv(n)-INyv(n),  Nyv+1(n),  .  .  Nyv+U1(n)  J, 

NZv^n)“l^zv(n^ *  Nzv+l(n)»  •  •  •»  Nzv+L-l(n)  1* 

Nxi(n),  Nyj(n)  and  Nzj(n)  are  introduced  here  for  simplicity  of  notation. 
They  are  given  by 

NxjCn)  -  N^n), 

Nyi(n)  -  DFT{ni(k-(i-l)T)}f 

and 

Nzi(n)  -  DFTtniOt-a-l)!)}. 

Six  matrices  Mj,  Nj ,  Pj,  Qj_,  Rj  and  S^  are  formed  vhere 


'  t 

1 

r 

1 

'  T  t 

|  | 

Mx(n)  - 

^(n)  .  . 

1 

•  *(m-L)<n> 

;  N^n)  - 

X2(n)  .  .  .  X(m_L+1)(n) 

1 

.  1 

1 

4 V 

1  1 

ir  1 
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< - 

4 - 

ft  t 

1  1 

Pl(n)  - 

Xl(n)  •  •  •  X(m-L)^n^ 

i  Qi(n)  - 

Y2(n)  .  .  .  Y(m_L+1)(n) 

1  ! 

.  i  i  j 

1  I 

i  i 

■  t  t 

1  1 

‘  T  t 

I  1 

Rj(n)  - 

Zj(n)  .  .  .  Z(m_L)(n) 

;  Sj(n)  - 

Z2(n)  .  .  .  Z(m_ul)(n) 

1  1 

.  1  i 

1  1 

.  1  l 

These  matrices  have  the  following  decompositions 

M^iO-Al  H(n)  S(n)  Ul(n)+Nl'(n)  ;  N^nl-Al  H(n)  ♦  Ul(n)+Nl"(n) ,  (4.2-22) 

P1(u)«A2  H(n)  S'  (n)  U2(n)+N2'  (n)  ;  Q^fO-AZ  B(n)  T  U2(n)+N2"(n) ,  (4.2-23) 

Rj(n)«A3  H(n)  S"(n)  U3(n)+N3' (n)  ;  B(n)  T  U3(n)+N3"(n) ,  (4.2-24) 

where  Ul,  U2,  U3,  S(n),  S'(n),  S'(n),  Nl'(n),  N2'(n),  N3'(n),  Nl"(n), 

N2"(n)  and  N3"(n)  are  given  by 

uiT-[  211,1  211,2  •  •  •  2Id,n<d>  1» 

U2T«[  221,1  2*1,2  *  •  •  22d,M<d> 

U3t-[  U3lfl  U3lf2  .  .  .  U3d,M<d>  ], 

21i,j  “I  1  ♦ij  -  *  (♦ij)(m_L+1)  ]T, 

U2ifj  .  [  1  Yij  •  •  (rij)<m-L+1>  ]T, 

21i,j  -  l  1  ♦ij  •  •  (♦ij)<®-L+1>  ]T, 

S(n)  -  diag{  S11(n)  S12(n)  .  .  .  SdM(d)(n)  }, 

S'  (n)  -  diag{  Su'(n)  S12'(n)  .  .  .  SdM(d)'(n)  ), 

S"(n)  -  diag{  Su"(n)  S12"(n)  .  .  .  SdM(d)"(n)  }, 

Nl'(n)-[NX1(n)  NX2(n)  .  .  .  NX(m_L)(n)  1, 

N2'(n)-[NY1(n)  NY2(n)  .  .  .  NY(m_L)(n)  ], 

N3'(n)-lNZ1(n)  NZ2(n)  .  .  .  NZ(m_L)(n)  ], 

Nl"(n)-[NX2(n)  NX3(n)  .  .  .  NX(m_ul)(n)  ] , 
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N2"(n)»[NY2(n)  NY3<n>  •  •  •  NY(m-L+l)<n) 

N3"(n)-[NZ2(n)  NZ3(n)  .  .  .  NZ(m_ul)(n)  ]. 

Assuming  the  signals  and  noise  to  be  statistically  independent  and  that  the 
noise  components  are  uncorrelated  from  sensor  to  sensor,  ve  get 

N-l 

E[M1H(n)M1(n)]  -  (1/N)  l  M1H(n)M1(n)  (4.2-25) 

n*0 

E[M1H(n)M1(n)]-UlH  VI  U1  ♦  L  N  a2  I(m_L)  , 

E[N1H(n)M1(n)]-Ula  t8  VI  U1  ♦  L  N  a2  Ii(m_L)  , 

E[P1B(n;P1(n)]«U2H  V2  U2  +  L  N  cr2  I(m_L)  , 

E[P1H(n)Q1(n)].U2H  I*  V2  U2  *  L  N  a2  I1(b_L)  , 

ElRj^nJR^n)]-^8  V3  U3  ♦  L  N  a2  I(m_L)  , 

E[S1H(n)R1(n)l-U3H  f8  V3  U3  +  L  N  a2  Ii(n_L)  , 
where  I(m_L)  is  the  (m-L)x(m-L)  identity  matrix  and  Ii(m_L)*  VI,  V2  and  V3 
are  the  matrices 

"0100.  .  .  .0 
0  0  1  0  ....  0 

0  0  0  1.  .  .  .0 

•  •  •  •  • 

Il(m-L)  •  •  •  •  • 

•  •  •  •  • 

0  0  0  0  ....  1 

.0000.  .  .  .0 

M1(n)-Al  H(n)  S(n)  Ul(n)+Nl'(n) 

Vl-E[SB(n)H8(n))A1H(n)A1(n)H(n)S(n)), 

V2«E(S'H(n)HH(n) )A18(n)A1(n)H(n)S' (n)] , 

Vl»Ef S"8(n)HH(n))A18(n)A1(n)H(n)S"(n) ] . 

Now,  let  M,  N,  P,  Q,  R  and  S  be  the  matrices 

M  -  ElM^Mj]  -  la2  I(b_L)  -  Ul8  VI  U1  , 

N  -  ElN^MiJ  -  La2  I1(o_L)  -Ul8  VI  Ul  , 
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(4.2-26) 

(4.2-27) 


(A. 2-28) 


P  -  E [ P !®P ! J  -  ha2  I(b_l)  -  U2a  V2  U2  , 

Q  -  ElQ^Pi]  -  Lcr2  Ii(ffl_L)  -U2H  T8  72  U2  ,  (4.2-29) 

R  -  E[R1HR1 ]  -  La2  I(m_L)  -  U3H  V3  U3  ,  (4.2-30) 

S  -  EJS^RjJ  -  La2  I1(m_L)  -U3h  f8  73  U3  .  (4.2-31) 

Consider  the  following  three  pencil  matrices  (M-AN),  (P-t|Q)  and  (R-vS). 

Note  that 

(M-XN).(Ul8  71  U1)-X(U1H  ^  71  U1)=U1H(I-X*H)71  U1  ,  (4.2-32) 

(P-HQ)«(U28  72  U2)-n(U2H  T8  72  U2)-U2HrT_r1rs)V2  U2  ,  (4.2-33 

(R-vS)«(U3h  73  U3)-v(U3H  t8  73  U3)*U3H(I-vta)73  U3  .  (4.2-34) 


The  matrices  Ul,  U2f  U3  i,  T,  and  T  are  all  of  rank  M  as  long  as  all  the 
^ij's,  the  Yij's,  the  +ij's  are  distinct  and  L>M.  Defining 

Epq,rs  =  ^  pq'^rs^* 

^  i,l 

L 

Epq,rs  *  ^  (y  pq"Yrs)^ 

F  i-1 

L 

Gpq,rs  *  pq-^rs^* 

clpq  ,  rs  *  ElS^n 

)Srs  (n)lH*^n)Hrs(n), 

C2pq,rs  -  E(S'^n)S'rs(n)  ]H${n)Hrs(n) , 
c^pq,rs  *  E[S"p^n)S"rs(n)  ]Hp^n)Hrs(n) t 
the  matrices  71,  72  and  73  can  be  written  as 


71- 


clllfllEll,ll 

c112,11e12,11 


L  ClM(d)M(d)fllEM(d)M(d)(11 


clll,M<d)M(d>Ell,M<d>M<d> 

cl12,n(d>M<d)E12,H(dVd> 


ClM(d)M(d)fM(d)M(d)EM(d)M(d)fM(d)M(d)  j 
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C211,11F11, 11 
c212,llF12,ll 


L  C2H(d)M(d)fllFM(d)M(d)fll 


C2llf„(d)M(d)FiifM(d)M<d) 

C2i2,M(d)M(d)F12>M(d)M(d) 


C2M(d)M(d)fM(d)M(d)FM(d)M(d)>M(d)H(d)  J 


c3ll,llGll,ll 

c312,llG12,ll 


[  C3M(d)M(d)fllGM(d)M(d)f n 


ClufM(d)M(d)Gii>M(d)H(d) 

C312,M(d)M(d)Gi2,M(d)M(d) 


C3M(d)M(d)>M(d)M(d)GM(d)M(d)tM(d)M(d)  J 


It  can  be  seen  that  the  matrices  VI,  V2  and  V3  are  of  rank  M.  The  rank  of 
the  pencil  (M-XN)  is  decreased  by  1  vhenever 

Xij-l/(4ij*)-(pij*)(A/cTs)sin(9j),  (4.2.35) 

for  i-1,  2,  .  .  .,  M<1>, 
j  > k-1 ,  2 ,  .  .  . ,  d . 

The  rank  of  the  pencil  (P-hQ)  is  decreased  by  1  whenever 

hij«l/(Yij*)-(Pij*)T,  (4.2.36) 

for  i-1,  2,  .  .  . ,  , 

j  ,  k—  1 ,  2 ,  .  .  . ,  d . 

The  rank  of  the  pencil  (R-vS)  is  decreased  by  1  whenever 

-ij-l/(*iJ*)-(PiJ*)<A/cTs)sin(eJ)  (Pij*)T>  (4.2-37) 

for  i-1,  2,  .  .  . ,  M^1), 
j , k— 1,  2,  •  .  • ,  d. 

As  was  noted  in  section  4.1,  from  the  first  set  of  generalized  eigenvalues, 


we  obtain  a  set  of  coupled  values  of  poles  and  angles  of  arrival.  The  sec¬ 
ond  set  solves  for  the  poles.  The  third  set  allovs  us  to  identify  which 
poles  go  with  which  angles  of  arrival  since 


vij  “  ^ij  •  *U j •  (4.2-38) 

This  way,  we  have  successfully  solved  for  the  angles  of  arrival  and  the 
poles  of  the  sources  using  a  model  of  wide  sense  stationary  signals. 
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4.2.1  COMPUTER  SIMULATION 

The  scenario  used  for  this  simulation  consisted  of  12  sensors 
uniformly  spaced  at  a  distance  A.  Again,  the  2  sources  were  assumed  to  be 
located  at  angles  0^»168  and  02*24°.  Two  cases  were  studied. 

Case  A 

In  this  case,  the  emitter  signals  are  assumed  to  have  been  genera¬ 
ted  by  passing  sequences  of  white  noise  through  linear  systems  with  impulse 
response  given  by  hl(k)-(pji)*c+(p^2)^  for  one  source  where  pn*0. 12+jO.  79 
and  Pi2*Pil*>  a°d  h2(k)-(p2j)k+(p22)k  for  the  other  source  with 
P2i*0.25+j0.90  and  P22“P21*’  A8a*n*  the  received  data  was  first  Fourier 
decomposed  using  128  snapshots  and  the  statistics  were  derived  using  50 
runs. 

Case  B 

In  this  case,  the  two  sources  are  assumed  to  have  identical  spec¬ 
tra.  The  emitter  signals  are  gernerated  by  passing  two  independent  white 
Gaussian  noise  sequences  through  a  linear  system  whose  impulse  response  is 
given  by 

h(k)-(p11)k+(p12)k 

where  Pn*0. 12+jO. 79  and  Pi2“Pll**  The  received  data  was  first  Fourier 
decomposed  using  128  snapshots  and  the  algorithm  vas  runs  50  times. 

The  results  are  plotted  in  figures  4.8  to  4.14  and  tables  4.7  to 
4.12.  In  these  plots,  (1)  denotes  the  estimates  for  Case  A  whereas  (2)  cor¬ 
responds  to  Case  B.  Both  methods  identify  the  poles  of  the  sources  with 
their  angles  of  arrival,  however,  riie  second  method  gives  better  esimates 
due  to  the  fact  that  the  poles  are  estimated  more  efficiently.  Ve  have  used 
the  smallest  window  possible  which  is  L-2  is  this  case  wheras  a  window  of 
length  L-4  was  used  in  the  first  case. 
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SAMPLE  MEAN 


SAMPLE  VARIANCE  (dB) 


SAMPLE  MEAN 


5 


ANGLE=24° 


(1) -No  Common  Poles 

(2) -  2Common  Poles 


Fig.  4.13.  Sample  Variance  of  Angle  at  24° 
Wide  Sense  Stationary  Signals. 
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MSE  (dB) 


ANCI-C=24  0 


(1) -No  Common  Poles 

(2) -  2  Common  Poles 


Fig.  4.14.  Mean-Squared  Error  of  Angle  at 
Wide  Sense  Stationary  Signals. 
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SNR  (dB) 

]  Sample  Mean 

I 

Sample  Variance 

30 

|  0.2499700-j0. 8999654 

1 

6. 5221090e-4 

25 

|  0. 2499475-j0. 8997493 

1 

1 . 1575040e-3 

20 

|  0. 2498654-j0. 8989769 

1 

2 . 0493467e-3 

15 

|  0. 2496969-j0. 8965538 

I 

3 . 6088992e-3 

10 

|  0. 2471662-j0. 8881544 

1 

1 

1 .0184543e-2 

5 

|  0. 1563175-j0. 8014570 

1 

5. 5371519e-2 

Table  4.7 

Sample  Mean  and  Variance  of  the  pole  P22’0. 25-jO. 90. 

(No  Common  poles) 


SNR  (dB) 

|  Sample  Mean 

1 

Sample  Variance 

30 

|  0 . 2500927+ jO . 8998703 

1 

5 . 1895110e-4 

25 

J  0. 250 1507 ♦jO. 8995793 

1 

1.0371592e-3 

20 

|  0. 250229 9+ j0. 8986719 

I 

1 . 8526319e-3 

15 

|  0. 2503528+j0 . 8959982 

1 

3.3094636e-3 

10 

|  0.2429177+jO. 8831745 

1 

1 . 9341080e-2 

5 

|  0. 1579504+j 0.8077933 

I 

5 . 6961089e-2 

Table  4.8 

Sample  Mean  and  Variance  of  the  pole  P21=0. 25*j0. 90. 

(No  Common  poles) 


1 18 


SNR  (dB) 

i  Sample  Mean 

1 

Sample  Variance 

30 

|  0.1200390-j0. 7898036 

1 

4 . 9322803e-4 

25 

|  0.1200731-j0. 7895529 

1 

8 . 7595923e-4 

20 

j  0.1201505-j0. 7889045 

1 . 5568603e-3 

15 

|  0.1204159-j0. 7871953 

! 

2 . 7756214e-3 

10 

|  0. 12 18340- jO. ^832218 

I 

4. 9832738e-3 

5 

|  0. 1287750-j0. 7795519 

1 

8 . 9591751e-3 

Table  4.9 

Sample  Mean  and  Variance  of  the  pole  p^2=0.  12- jO.  79. 

(No  Common  poles) 


SNR  (dB) 

1 

Sample  Mean 

1 

Sample  Variance 

30 

1 

0. 1199508+jO. 7900655 

1 

4 . 9638440e-4 

_ 

A 

1 

0.1199184+j0. 7900197 

1 

8 . 8051870e-4 

20 

1 

0. 1 1 98816+ jO . 7897383 

1 

1 . 5633145e-3 

15 

1 

0. 1199524+jO.  7886938 

1 

2 . 7888690e-3 

10 

1 

0. 1210258+j0. 7859387 

1 

5.0294539e-3 

5 

1 

0.1271952+j 0.7843727 

1 

8 . 9402702e-3 

Table  4.10 

Sample  Mean  and  Variance  of  the  pole  Pn=0. 12+j0. 79. 

(No  Common  poles) 
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SNR  (dB) 

j  Sample  Mean 

1 

Sample  Variance 

30 

|  0.1200018-j0. 7899857 

1 

1.0762867e-4 

25 

|  0.1200033-j0. 7899747 

1 . 9138129e-4 

20 

|  0.1200061-j0. 7899654 

I 

3 . 4014991e-4 

15 

|  0.1200114-j0. 7899263 

I 

6 . 0441019e-4 

10 

|  0.1200223-j0. 7898808 

1 

1.0735721e-3 

5 

|  0.1200452-j0. 7898291 

1 

1 . 9071890e-3 

Table  4.11 

Sample  Mean  and  Variance  of  the  pole  Pi2=P22=® ■ 12-jO • 79 

(2  Common  poles) 


SNR  (dB) 

|  Sample  Mean 

1 

Sample  Variance 

30 

|  0. 1199916+j0. 7900093 

1 

1 . 3058618e-4 

25 

|  0.1199852+j0. 7900167 

1 

2 . 3202586e-4 

20 

|  0.1199739+j0. 7900308 

i 

4. 1216239e-4 

15 

|  0.1199542+j0. 7900584 

i 

7 . 3149381e-4 

10 

|  0.1199203+j0. 7901148 

1 

1 . 2967319e-3 

5 

|  0.1198631-t-jO.  7902434 

1 

2 . 2957609e-3 

Table  4.12 

Sample  Mean  and  Variance  of  the  pole  Pn=P21-0-  12+jO.  79. 

(2  Common  poles) 
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4.3  FOURIER  APPROACH 


Consider  the  problem  of  estimating  the  angles  of  arrival  of 
videband  signals.  The  notion  of  Fourier  coefficients  is  used  here  in  con¬ 
junction  with  the  matrix  pencil  approach.  Assume  that  all  incoming  signals 
have  approximately  the  same  bandwidth  B.  Let  T  be  an  observation  interval 
and  denote  by  and  tog,  the  lowest  and  highest  frequencies  contained  in  B. 
In  practice  the  frequency  content  is  determined  by  Fourier  decomposition  of 
the  received  signals.  The  received  signal  at  the  i-th  sensor  can  be  modeled 
as 

d 

x1(t)-E  8(0^)  sk(t-Tik)  +nj(t);  i-1,2,.  .  .,m,  (4. 3-i) 

k-1 


where  t ^  is 


Tik=*  <  *  ~ 1 )(  A/  c)  s  in(  8^) . 
Define  the  Fourier  coefficients  as 


*!<«*)- 


(T)* 


T 


/2 

XiCOexpf-jWrt}  dt  , 
T/2 


(4.3-2) 


where  R  is  the  number  of  subbands,  A6>-(ug>u^)/R>iWidth  of  each  subband, 
«r»(2it/T)(r^+r),  rj  is  a  suitably  chosen  integer  such  that 
(2ii/T)r1-(<uL+/W2)  and  (2n/T)(r1+R)-(uB-A«/2) . 

Taking  the  Fourier  coefficients  of  both  sides  of  equation  (4.3-1), 

we  get 


d 

Xi(«r)  -  Z  a(ek)Sk(«r)e^i-1>*k(“r>  ♦  Ni(a»r);i«l,2, .  .  .,m,  (4.3-3) 

k-1  r*l,2,.  .  .,R 


where 
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(4.3-4) 


♦jt(wr)«-(ttr)(  A/c)sin(  0^) . 

Given  this  set  of  Fourier  coefficients,  (m-L+1)  vectors  of  length  L 

are  formed  vhere 

d  <  L  <  (m-d), 

X^MW  Xn+1(*t)  •  •  ♦  Xn+L_1(«r)JT  ;  n-1,  2 . (m-Ul) . 

r*l,  2,  .  .  . ,  R 

It  can  be  shown  that  Xj^Qj.)  can  be  put  in  the  form 

Xn(«r)-A(«r)*(n-1>(<dr)BS(a>r)+Nn(a)r),  (4.3-5) 

vhere 


A(«t> 


i)*l<“r> 

• 

• 

•  • 

•  • 

.  1 

.  ej*d(“r) 

•  • 

•  • 

♦ 

^(L-D^Cov) 

•  • 

^)(L-l)#2(»r)  ! 

•  • 

;  ej(L-Dfd(tor) 

*  a(  0[^) , 

B  -  diag  {a^  &2  •  •  -a<i  }» 

♦(«,.)  »  diag  {  ej*l<“r>  eJ+2<*r>  .  .  .ej*d<“r>  ), 

SCttj.)  -  I S2<*>r)  S2(«r)  •  •  •  S,j(«r)  1 

and 

Hn^(®r)  *  1^(0^)  •  *  *  ^n+L-l^^r^l* 

Non  singular  transformation  matrices  Tj^Wp)  of  dimension  (LxL)  are  then 

used  in  such  a  way  that 

Tn(o»r)A(«r  )#(n-1><«or>-A(«u0)  *(n~1>(at>) ,  (4.3-6) 

vhere  is  a  conveniently  chosen  frequency.  In  this  fashion  all  the  power 
in  the  corresponding  sub  bands  is  "moved"  to  a  single  band.  This  process  is 
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lenaed  in  the  literature  as  focusing.  It  is  desirable  to  solve  equation 
(4.3-6)  for  T^toj.).  However,  («£.)  is  of  dimension  (Lxd)  and, 

therefore,  does  not  possess  an  inverse.  Vithout  loss  of  generality,  it  is 
possible  to  augment  A(a>r)*(n_^(a)r)  by  a  matrix  Vfaij.)  of  dimension 
(Lx(L-d))  so  as  to  generate  the  non  singular  square  (LxL)  matrix 
(A(«r)*(n-1)(wr)  «(«,.)]. 

At  c^),  this  matrix  becomes 

[A(aJ0)*Cn-1>(030)  W(«o)  ] . 

An  equivalent  equation  for  equation  (4.3-6)  is  then  given  by 
Tn(Wr)[A(«r)*(n-1)(«r)  V(Wr)]-[A(u0)*(n-1)(U0)  U(UQ)]. 

It  follows  that 

Tn(<or)-[A(wD)*(n-1)(t^)  W(^)][A(«r)*<n-1)(o>r)  W(«r)]~1. 

Let  &i,  02*  *  *  •»  3d  be  preliminary  estimates  of  the  angles  of  arrival  ©j, 
©2»  •  •  •*  3d  obtained  by  some  simple  low  resolution  technique  such  as  the 
periodogram.  Define 

♦k  («!■)— («r)(A/c)sin(0jt). 

AtWj.)  is  then  approximated  by  the  matrix 


f  exp{j41(a»r)} 


A3(“r)»  • 


.  1 

exp{j4d(“r)) 


[  exp{j(L-l)f1(»r)} 


exp{j(L-l)4d(o>r)) 


Ap(o^)  is  obtained  from  A^u^)  by  replacing  Wj.  with  The  desired  trans¬ 
formations  are  then  given  approximately  by 

Tn(«r)-lA0(«o)*0(n"1)(wO)  w(«D)](A0(wr)*3(n‘1)(a>r)  ^Vl"1  ■ 

To  prevent  matrices  from  being  singular,  V  is  chosen  to  have  the  same  form 
as  A  but  is  evaluated  at  distinct  angles  different  from  Hj,  32»  •  •  •»  3d’ 
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If  it  happens  that  all  the  true  angles  of  arrival  are  within  the  neighbor¬ 
hood  of  a  single  angle  0,  the  approximate  transformation  matrices  are 
diagonal  and  are  of  the  form 

Tn(«r)-e-j(n-1X“0-«rHA/c>sin<®)  T^) 

where 

T1(»r)-diag{  1  «-J<**O-«t>(^Osin(0)  .  .  .  «-j(I*--X«*)-“r><A/c)sin(0)j  # 

Applying  these  transformation  to  every  vector,  we  obtain 

Tn(«r)Xn(«r)  -  Tn(o»r)A(«r)*<n-1>(«r)BS(«r)  +  Tn(«r)Nn(«r) 

»  A(o^)*(n-1)(u0)BS(«r)  +  Tn(«r)Nn(o)r).  (4.3-7) 
With  respect  to  the  R  sub-bands,  consider 

R 

*„(<■*))  -  (1/R)  l  Tn(oar)XnC«r). 

r-1 

Let  S'  and  be 

R 

S'  -  (1/R)  l  S( tOj- )  =  [S' !  S'2  •  •  •  S'd  JT, 
r*l 

R 

Bn  -  (1/R)  £  V“r)Bn<“r>-i"  n  "'n.L-1  1T- 

r-1 

Therefore,  Xj^o^))  can  be  expressed  as 

Xn(«0)-A((^)*<n-1)(«^))BS'+N^  .  (4.3-11) 

In  the  remainder  of  this  discussion  the  dependence  on  is  assumed.  The 
two  matrices  and  N}  are  then  formed  where 


'  T  T  T 

1  1  1 

4 - 

4 - 

4 - 

Mi  - 

l\l2  •••  X(m-L) 

;  Nj  - 

ll  *3  • • '^(m-L+1) 

4 - 

4 - 

« - 

1 _ 

1  1  1 

l  <1  4. 

(4.3-8) 


(4.3-9) 


(4.3-10) 
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These  can  be  decomposed  as 


t  t 


Mx  «  ABS'  AB*S'  .  .  .  AB*(m-L_1)S'  +  N' 


(4.3-12) 


4.  4, 


T  t 


Nj  «  AB*S'  AB^S'  .  .  .  AB*(®-L^S'  +  N", 


(4.3-13) 


4  l 


where 


T  T 

I  I 


N'  -  Nj'  .  •  .  N(a.L)' 


I  I 

l  1 
T  T 


N2'  «3' 

I  I 

1  4. 


N(m-Ul) 


Simplification  of  and  results  in 

Mj  -  AB[ S'  *(®-L-l)s']  +  N', 

Nx  .  AB*[S'  ^S'  .  .  :  *(®-L-l)s']  +  N" 


Let  F  be  the  matrix 


F  -  S  *S  .  .  .  *(®-L-l)s 
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The  matrix  P  can  be  written  as 
P.DU, 

where 

D  •  diag{  S'i  S'2  •  •  •  S'd  }, 

r  i  ei*i  .  .  .  eK®-L-mi 
1  eJ+2  .  .  .  ej(m-L-1)^2 

U  -  .  . 

.  1  eJ*d  .  .  .  ej(m"L-1)+d 

Then 

-  ABDU+N'  (4.3-14) 

and 

Nx  »  ABD*U+n".  (4.3-15) 

Assuming  that  the  signals  and  noise  are  statistically  independent  and  that 
the  noise  components  are  uncorrelated  from  sensor  to  sensor,  we  get 

J-UHVU  +  E[N'HN' ]  (4.3-16) 

E(N1HM1]-UH*HVU  ♦  E[N"HN']  (4.3-17) 

where  V  is  the  matrix  V-E[DHBHAHABD1 .  Defining 

Fpa  »  E  ( i-1) ( ♦p-fq) f 
i-1 

spq  -  ElSqSp  ], 
apq  -  aJsp» 

the  matrix  V  becomes 

'  sllallpll  ....  SdladlPdl 

s12a12p12  ....  Sd2ad2Pd2 

V  - 

■  sldaldFld  ..••  SddaddFdd  .  . 

Note  that  the  matrix  V  is  of  rank  d  even  in  the  presence  of  fully  corre¬ 
lated  sources.  Define  the  matrices  N  and  N  as 

M  .  EIM1HN1]-EIN,HN' ]  -  UH  V  U  (4.3-18) 
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and 


N  -  E[N1hM1]-E[N"hN'  )  .  UH  ♦“  V  U.  (4.3-19) 

The  matrix  pencil  then  becomes 

M-XN-UHVU~XUH*HVU-UH(I-X#H)VU  (4.3-20) 

which  satisfies  the  requirements  of  the  pencil  theorem.  Hence,  the  values 
of  X  for  which  the  rank  of  M-XN  decreases  by  1  are  given  by 

\  -  eJ*k  ;  k-1,2 . d.  (4.3-21) 

The  angles  of  arrival  are  given  by 

8^  -  sin"1(jcln(Xjt)/ci^h) ;  k-1,2,  «..,d.  (4.3-22) 


4.3.1  SIMULATION  3 

Several  possibilities  exist  for  choosing  the  transformation 
matrices  Tn  [581.  It  can  be  shown  that  a  diagonal  transformation  leads  to 
the  simplest  analysis.  Assuming  the  sources  to  be  clustered  within  the 
proximity  of  one  location  0,  the  transformation  matrices  Tn ( >  then  become 

Tx ( Mj. ) 

where 

T^foic )-diag{  1  e~J  ^_,ar A/c)sin( (5)  .  .  .  e-j  (L-l )  (o^-ca,. )( A/c  )sin(  (3)  j 
With  this  transformation  it  follows  that 

Tn(«r)A(o»r)#(n-1)(«r)  -  A(^)#(n-D(a^). 

Assuming  that  the  noise  components  are  uncorrelated  from  sensor  to  sensor 
and  from  sub-band  to  sub-band  with  zero  mean  and  variance  <y*  ,  it  can  be 
shown  that 

E[N'hN' 1  -  R  a2  I(m.L)d 
E[N"«N'l  ♦  R  a2  Il(m_L) 

where  I(m_L)  is  the  (m-L)x(m-L)  identity  matrix  and  Il^m_L)  f^e  matrix 


0  1  0  0  ....  0  • 

0  0  1  0  ....  0 

0  0  0  1  ....  0 


•  «  •  «  • 

0  0  0  0  ....  1 
.0000.  .  .  .oj. 

In  the  simulations,  ve  h.'.ve  considered  a  linear  array  consisting 
of  8  sensors  uniformly  spaced  at  a  distance  &-c/(2  fg)  .  Following  the  ex¬ 
ample  in  [57],  the  two  sources  were  assumed  to  be  located  at  16°  and  24° 
and  to  have  ideal  rectangular  spectra  of  bandwidth  B-40  Hz  centered  at 
f0«100  Hz.  The  broadband  signals  were  first  decomposed  into  33  narrowband 
components.  100  snapshots  were  taken  for  each  of  the  50  runs-  As  in  chapter 
3,  ESPRIT  can  be  used  either  with  overlapping  subarrays  or  non  overlapping 
subarrays.  In  the  first  case,  the  subarray  X  consists  of  the  first  7 
sensors  and  the  subarray  Y  consists  of  the  last  7  sensors.  In  the  second 
case,  two  adjacent  sensors  were  considered  as  a  pair.  The  results  of  the 
simulation  are  plotted  in  figures  4.15  to  4.20.  In  these  figures,  the 
Moving  Window  is  represented  by  (1),  ESPRIT  overlapping  case  by  (2)  and 
ESPRIT  non  overlapping  case  by  (3).  (4)  represents  the  Cramer-Rao  lower 
bound  (CRLB)  vhich  is  described  in  the  appendix. 

Note  that  the  estimates  obtained  through  moving  vindow  have  small 
bias  and  their  variances  approach  the  CRLB  very  closely  especially  at  high 
SNR.  Ve  have  thus  shown  that  the  moving  window  can  be  applied  in  conjunc¬ 
tion  vith  CSS  and  that  it  performs  slightly  better  than  ESPRIT. 
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SAMPLE  MEAN 


SNR(dB) 

( 1 )  -Mov ing  Window 

(2) -ESPRIT:  Linear  Overlapping  Case 

(3) -ESPRIT:  Linear  Non  Overlapping  Case 


Fig.  4.15.  Sample  Mean  of  Angle  at  16 
Fourier  Approach. 


1 29 


-10  -5  0  5  10  15  20  25 

SNR(dB) 


(1) -Moving  Window 

(2) -ESPRIT:  Linear  Overlapping  Case 

(3) -ESPRIT:  Linear  Non  Overlapping  Case 

(4) -CRLB 

Fig.  4.16.  Sample  Variance  of  Angle  Estimate  at  16 
Fourier  Approach. 
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MEAN-SQUARED  ERROR  (dB) 


SNR  (dl3) 


(1) -Moving  Window 

(2) -ESPRIT:  Linear  Overlapping  Case 

(3) -ESPRIT:  Linear  Non  Overlapping  Case 

(4) -CRLB 

Fig.  4.17.  Mean-Squared  Error  of  Angle  at  16 
Fourier  Approach. 
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SAMPLE  MEAN 


ANGLE=24° 


(1)  -Moving  Window 

(2) -ESPRIT:  Linear  Overlapping  Case 

(3) -ESPRIT:  Linear  Non  Overlapping  Case 


Fig.  A. 18.  Sample  Mean  of  Angle  Estimate  at  24° 
Fourier  Approach. 
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(1) -Moving  Window 

(2)  -ESPRIT:  Linear  Overlapping  Case 

(3) -ESPRIT:  Linear  Non  Overlapping  Case 

(4) -CRLB 

Fig.  4.19.  Sample  Variance  of  Angle  at  24° 
Fourier  Approach. 


MEAN-SQUARED  ERROR  (dB) 


ANGLE=24° 


(1) -Moving  Window 

(2)  -ESPRIT:  Linear  Overlapping  Case 

(3)  -ESPRIT:  Linear  Non  Overlapping  Case 

(4) -CRLB 

Fig.  4.20.  Mean-Squared  Error  of  Angle  at  24° 
Fourier  Approach. 
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CHAPTER  5 


PERTURBATION  ANALYSIS 


The  methods  described  previously  assume  that  additive  noise  can  be 
suppressed  through  noise  compensation.  Also,  the  uniform  spacing  between 
the  sensor  eleaents  of  an  array  is  assumed  to  be  knovn.  In  practice,  how¬ 
ever,  it  is  likely  that  the  noise  compensation  will  be  non  ideal  and  that 
the  sensor  elements  will  be  perturbed  from  their  uniform  spacing.  In  this 
section,  performance  degradation  is  investigated  due  to  imperfect  compensa¬ 
tion  of  the  additive  noise.  The  case  of  offsets  in  the  sensor  spacing  is 
studied  in  a  similar  fashion.  The  chordal  metric  [85]  is  introduced  as  a 
measure  of  the  distance  between  the  true  and  perturbed  eigenvalues. 
Theoretical  upper  bounds  are  derived  for  the  chordal  metric  for  both  the 
Moving  Vindov  and  ESPRIT. 

5.1  Chordal  Metric 


Let  C  denote  the  field  of  all  complex  numbers.  Consider  the  eigen¬ 
value  problems 


M  x  -  X  N  x 


(5.1-1) 


and 

M  -  A  £H  N  (5.1-2) 

vhere  H  denotes  complex  conjugate  transpose,  x  and  jr  are  called  the  right 
and  left  eigenvector,  respectively,  of  the  pencil  formed  by  M  and  N.  Solu¬ 
tion  for  £  proceeds  by  solving 

Mh  x  -  A*  Nh  y  . 

Introduce  the  Euclidian  matrix  norm  defined  as 


|  |M  j  j  -  sup  (  jMx  I  |  -  sup{x®(M®M)x}'*. 
Mil  l-i  ll*ll-i 
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Ve  are  interested  in  the  generalized  eigenvalue  problem 


where 


M  x  ■  X  N  x 

M-M  +  AM-M  +  E 
N-N+AN-N  +  F. 


(5.1-3) 


(5.1-4) 

(5.1-5) 


Let  04  and  be  the  quantity 

04  -  (5.1-6) 

and 

•  ZiHf,5i*  (5.1-7) 

where  xj  and  £4  are  the  i-th  right  and  the  i-th  left  eigenvectors  of  the 
pencil  formed  by  M  and  N.  it  follows  from  equation  (5.1-1)  that 


Xi  «  04/3^.  (5  1-8) 

Stevart  [85]  showed  that  snail  perturbations  in  M  and  N  result  in 


a4+^i3Exj+0(e^)  c4'+0(e2) 

Xj  .  - - - —  ■  - r~  (5.1-9) 

3i+Xi0F*i+O(c2)  Pj'^e2) 

where 

0(e2) 

lin  -  -  0. 

e-O  c 

Define  the  chordal  aetric  as 

I  h-h  I 

XCXj,^)  .  -  .  (5.1-10) 

J  j  i.|XjI2' 

With  this  definition  it  was  shown  [85]  that 

X(Xi,Xi)  <  c/Yi  +  (He2)  (5.1-11) 

where 

t  -J  ||Ej|2  ♦  |  |F  |  |2  (5.1-12) 

and 
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Ti  «J  <42  +  3i2 


(5.1-13) 


In  our  applications,  the  eigenvalue  Xj  is  related  to  the  angle  of 
arrival  0^  through  the  equation 

exp{j(wA/c)sin(@i)) . 

The  perturbed  eigenvalue  Xj  then  corresponds  to  an  angle  of  arrival  0^ 
given  by 

Xj»  exp{j (»A/c)sin(0i)} . 

Let  ♦j*(ttA/c)sin(0j)  and  ♦i»(«A/c)sin(0j) .  It  can  be  shown  that 

|  Xj-Xj  |  -  2  sinU^-**)^}.  (5.1-14) 

Note  that  | |Xj | |«1  and  ||X^||-1.  Therefore,  equation  (5.1.10)  reduces  to 

XfXj.Xi)  -  sinK^-^)^}.  (5.1-15) 

which  implies  that 

easin'1  {sinfe^  ±  (2c/wh)sin'1tx(Xi,Xi)  H  .  (5.1-16) 

Hence,  given  the  value  of  the  chordal  metric,  it  is  possible  :o  determine 
the  perturbed  angle  of  arrival  by  using  equation  (5.1-16). 

5.2  PKRTUBATION  DUE  TO  NOISE 

In  this  section  we  study  the  effects  of  non  ideal  noise  compensa¬ 
tion  on  the  performance  of  the  Matrix  Pencil  Approach.  To  compensate  for 
the  noise,  it  was  shown  in  chapter  2  that  it  is  necessary  to  know  the  noise 
covariance  matrix.  However,  in  pactice,  the  noise  covariance  matrix  is  not 
known  exactly.  To  obtain  a  measure  of  the  perturbation  in  the  eigenvalues, 
upper  bounds  are  derived  for  both  the  moving  vindov  and  ESPRIT  operators. 
5.2.1  MOVING  VINDOV 

For  the  Moving  Vindov  operator  discussed  in  chapter  2,  two 
matrices  E{MjHMj]  and  E(NjHM^]  are  formed  where 

B[M1hN1J-UhVU  ♦  Lo2  I(B_l)  (5.2. 1-1) 
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B[N1®M11»U®#®VU  +  La2  I1(a_L)  (5. 2. 1-2) 

wher«  U,  V  and  UH  are  (m-L)x(m-L)  non-singular  matrices  and  I(m_L)  is  the 
(m-L)x(m-L)  identity  matrix  and  Ii(m_L)  is  the  (m-L)x(m-l)  matrix 

'0100.  .  .  .  0  ' 

0  0  1  0  ....  0 

0  0  0  1.  .  .  .0 

•  •  •  •  • 

^l(m-L)  *  *  •  •  • 

•  ■  •  •  • 

0  0  0  0  ....  1 

.0000.  .  .  .  0 

Let 

M  -  EIM^] 

N  -  ElN^Mj] 

m«uhvu 
E-Lff2  I(B_l) 

N-UHf®VU 
F-Lo2  I1(b_L). 

As  a  vorst  case,  assume  noise  correction  is  not  attempted.  Thus, 

M-M+AM-M+E  (5. 2. 1-4) 

N-N  +  AN-N  +  F.  (5.2. 1.5) 

In  order  to  use  the  bound  on  the  chordal  metric,  one  needs  to  evaluate  the 
euclidean  norms  of  the  matrices  E  and  F.  for  this  purpose,  consider 
| |Ex| |2  -  {(Ex)H(Ex)}-  xh(EhE)x  . 

Maximizing  ||Bx||  is  the  same  as  maximizing  the  quadratic  form  xH(EHE)x. 
Given  a  Heraitian  matrix  M,  it  is  knovn  that  the  maximum  value  of  the  quo¬ 
tient, 

xH  M  x 


is  equal  to  the  largest  eigenvalue  of  M.  Therefore,  under  the  constraint 
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xHx-l,  the  maximum  value  of  xH(EHE)x  is  equal  to  the  maximum  eigenvalue  of 
E®B.  Since  B-  Lo2  I(a_L)  »  then  E^E  ■  L2  o*  I(m-L)’  The  largest  eigenvalue 
of  this  matrix  is  L2  a A.  Bence, 

|  |E  |  |-L  a2-  (5. 2. 1-6) 
Similarly,  because  F»Lff^  Il(m-L)*  F®F  is  t^ie  following  matrix 


pHp 


L2  a* 


r  o  o 
o  1 


o  o 
o  o 


o  o 


1  o 


L  o  o 


o  1  J 


The  largest  eigenvalue  of  this  matrix  is  also  L2  o**.  Thus 

||F||«L  a2.  (5.2. 1-7) 

Therefore, 


e  J  |  |E|  |2  +  |  |F|  |2  »L  a2  fT  .  (5.2. 1-8) 

For  e-*0  ;  i.e,  o2<<(l/L(2;^!) ,  the  bound  on  the  chordal  metric  then  becomes 


X(Xi,Xi)<  L  a2 


<<Zi  M*i)Z+(Zi  N5i>Z) 


(5. 2. 1-9) 


5.2.2  ESPRIT 

ESPRIT  can  be  employed  in  a  variety  of  situations.  The  general 
case  involves  isolated  doublets  located  randomly  in  the  plane.  However, 
when  a  linear  uniformly  spaced  array  is  used,  two  schemes  are  possible 
depending  upon  vhether  the  doublets  overlap  or  not.  For  a  given  number  of 
sensor  elements,  the  overlapping  case  has  the  advantage  of  having  a  larger 
aperture  size  and  thus  a  better  resolution. 
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5. 2. 2.1  General  Case 


For  the  general  case  of  ESPRIT  presented  in  chapter  2,  tvo 
matrices  were  formed  from  the  data  vectors  X  and  Y  such  that 

E[ X  X^-A^0  +  a2  Im  (5. 2. 2. 1-1) 

E[X  .  (5. 2. 2. 1-2) 

Let 

M  »  E[X  X0] 

N  «  E[X  Y0] 

M-AjSA!3 

E-^  Im 

N.AjS^Aj0 

F*0. 

Thus,  assuming  noise  correction  is  not  attempted, 

M-M+AM.H  +  E 
N.N+AN.N  +  F. 

Ve  know  that  ||E||  and  )  | F  |  |  are  equal  to  the  square  root  of  the  largest 
eigenvalue  of  E0E  and  F0F,  respectively.  Since  E«  a2  Im  and  F*0,  EHE=a4Im 
and  F0F»O.  The  largest  eigenvalue  of  these  matrices  are  and  0.  There¬ 
fore, 

11*11-  o2  (5.2.2. 1-3) 

||F||-0.  (5. 2. 2. 1-4) 

Thus, 

C  -J  1  |E  I  I2  7  I  |F  I  I2  .  a2  .  (5. 2. 2. 1-5) 

For  e  -»  0,  the  bound  on  the  chordal  metric  becomes 
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X(XifXi)^ 


(5.2.2. 1-6 


l((£iHNXi)2*(£iHNXi)2) 


5. 2. 2. 2  Linear  Array:  Overlapping  Doublets  Case 

In  this  case  a  linear  array  composed  of  m  sensors  is  used.  1'he 
signal  received  at  the  i-th  sensor  can  be  modeled  as 


yi<t,6)«  I  ajlSic(t)  e-^*  ♦  n«(t)  ;i-l,  2, 

k-1 


.  , ,m,  (5. 2. 2. 2-1) 


where 


A|t«(«A/c)sin(Ojt)  ,  k-1,  2, 


,2. 2. 2-2) 


nj(t)  is  the  additive  noise. 

Two  overlapping  subarrays  Yj  and  Y2  are  then  formed  where 

II  -  t  y'l  Y2  •  •  •  Y(m-l)  1T  <- 

Y2  -  t  y2  Y3  •  •  •  Ym  1T-  (' 

In  chapter  3,  it  was  shown  that 

K f XlIlH  1  -A2BSBhA2h  ♦  a2  (! 

E(YiY2Hl-A2BS*HBHA2H  *  ^  I2(m_i)  ( - 

where  I(m_i)  is  the  (m-l)x(m-l)  identity  matrix  and  I2(m_i)is  the 
(m-l)x(m-l)  matrix  shown  below 


.2. 2. 2-3) 
.2. 2. 2-4) 


2. 2. 2-5) 


2. 2-6) 


0  0  0  0 
10  0  0 
0  10  0 


rKm-u 


0  0  0  0 
0  0  0  0 
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M  .  ElY^l 
N  -  ElY^l 

m.a2bsbha2h 

E“ff2  *(m-l) 

n-a2bs^bha2h 

?ma2  Il(m-1)- 

Thus,  assuming  no  compensation  for  the  noise, 

M-N-t-AN'.M  +  E 
N  •  N  >  I  •  N  i  F. 

As  before,  we  know  that  ||E||  and  ||F||  are  equal  to  the  square  root  of  the 
largest  eigenvalue  of  EHE  and  FHF,  respectively.  It  follows  that 

!  IKI  I-  ff2  (5. 2. 2. 2-7) 

and 

I  lF  I  I-  02-  (5. 2.2. 2-8) 

Thus, 


t  |  |E|  |2  ♦  |  |F|  |2  -  17  a2  , 


(5. 2. 2. 2-9) 


For  c  +•  0,  the  bound  on  the  chordal  metric  becomes 


X(  Aj ,  Aj )  <|  e2 


((jriHMXi)2  +  (£iHNxi)2) 


(5.2.2.2-10) 


5. 2. 2. 3  Linear  Array:  Non-Overlapping  Doublets  Case 

In  this  case  two  non-overlapping  subarrays  Yj  and  Y2  are  generated 
from  the  data  received  at  the  sensor  array.  Assuming  m  to  be  even,  Y^  and 
Y2  are  given  by 

Xl  -  l  Vl  Y3  •  •  '  y(m-l)  1T  (5.2.2. 3-1) 


U2 


(5. 2. 2. 3-2) 


h  -  t  y2  y4  •  •  •  y«  IT- 

As  shovn  in  chapter  3 

BlMlH)-*3BSBlV  *  ^  hm/2) 

E(YiY2H]-aiBs*HbHaiH  . 

Let 

M  -  ElYiYiHl 
N  «  E[ Y1Y2h] 

M-a3BSBHa3H 

S-*2  Im 

N-a3BSBh*Ha2H 

F-O. 


Thus 

N-N  +  AM-N  +  E 
N-N+AN-N  +  F. 


|  |E  |  |  and  |  |F  |  |  are  evaluated  as  before.  Since  E-o2^  and 

F-0,  EHE.„*I(m/2) 

and  FhF-0.  The  largest  eigenvalue  of  these  matrices  are  a * 

and  0,  There- 

fore, 

||E||.  J 

(5. 2. 2. 3-5) 

and 

||F||.  0. 

(5. 2. 2. 3-6) 

Thus, 

t.r 

1  |E|  I2  *  1  |F|  I2  -  ^  ■ 

(5. 2. 2. 3-7) 

In  this  case,  the 

bound  on  the  chordal  metric  becomes 

X(Xj,Xj)^  ■ 

*2 

(5. 2. 2. 3-8) 

• 

1 

l((£iHMxi)2*(JriHNxi)2) 

14  3 

(5. 2. 2. 3-3) 
(5, 2. 2. 3. 4) 


5.2.3  COMPUTER  SIMULATION 

A  computer  simulation  was  carried  out  to  demonstrate  the  ap¬ 
plicability  of  the  upper  bounds  derived  for  the  chordal  metric.  It  should 
be  pointed  out  that  fev  adjustments  had  to  be  made  in  order  to  exactly  use 
the  derived  bounds.  The  bounds  involve  the  matrices  M  and  N  and  their  cor¬ 
responding  eigenvalues  and  eigenvectors,  the  latter  being  of  dimension 
((m-L)xl).  However,  from  previous  sections,  we  have  seen  that  the 
dimensionas  of  N  and  N  are  reduced  to  the  order  needed  (d,  the  number  of 
sources)  so  that  only  the  signal  eigenvalues  are  estimated.  The  IMSL 
routine  BIGZC  called  to  do  this  will  return  eigenvectors  of  dimension  dxl. 
There  exists  methods  to  obtain  the  desired  eigenvectors  from  the  returned 
set.  Ve  opt  for  the  following.  Recall  that  the  original  problem  involved 
solving  the  equation 

Mx  -  X  Nx.  (5.2.3-n 

The  singular  value  decomposition  (SVD)  of  the  matrix  N  results  in 

Let  N+  be  the  pseudo  inverse  of  N.  N+  satisfies  the  Moore  Penrose  equa¬ 
tions.  It  is  clear  that  N+  is  given  by 
N*.Vn  (S,,)-1  UnH, 

where  (Sn)-^  consists  on  the  inverse  of  the  non  zeros  singular  values.  Pre¬ 
multiplying  both  sides  of  equation  (5. 2. 3-1)  by  N+  results  in 

N+M  x  -  X  N+N  x.  (5. 2. 3-2) 

Noting  that  N*N  is  the  identity  matrix  provided  that  x  is  in  the  range  of 
N,  solution  of  equation  (5. 2. 3-1)  is  equivalent  to  solving  the  equation 

N+M  x  -  X  x.  (5. 2. 3-3) 

The  eigenvalues  and  eigenvectors  obtained  in  this  fashion  would  be  of 


dimension  (m-L)xl  . 

The  scenario  used  for  the  simulation  consisted  of  a  linear  array 
of  8  sensors  uniformly  spaced  at  a  distance  A.  2  incoherent  sources  are 
present  and  are  located  at  angles  9j«16°  and  .  The  incoherent  case 

was  chosen  so  as  to  give  a  fair  comparison  to  ESPRIT.  The  additive  noise 
vas  generated  as  white  Gaussian  with  zero  mean  and  unit  variance.  The 
perturbed  angles  of  arrival  of  the  sources  are  obtained  using  the  upper 
bound,  the  chordal  metric  and  the  perturbed  eigenvalue.  Tables  5.1  to  5.6 
shov  the  sample  mean  of  the  angle  estimates  obtained  from  50  runs  where  100 
snapshots  were  considered  in  each  run.  Note  from  these  tables  that  the 
bounds  derived  in  this  section  perform  quite  well  compared  to  the  exact 
value  of  the  chordal  metric.  In  some  instances,  the  bound  is  even  smaller. 
This  due  to  the  bias  that  exists  in  the  magnitude  of  the  eigenvalue. 
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Moving  Vindov 


Angle-24® 


SNR 

(dB) 

9  obtained 
from  Bound 

9  obtained  from 
chordal  metric 

8  obt§ined 
from  X 

10 

24.26962 

24.56922 

24.13065 

9 

24.45502 

24.65348 

24.14139 

8 

24.71086 

24.73565 

24.15029 

Table  1 
(Sample  Mean) 


Angle-160 


SNR 

(dB) 

9  obtained 
from  Bound 

9  obtained  from 
chordal  metric 

9  obtjined  1 
from  X  | 

10 

|  16.56222 

|  16.48074  | 

|  16.00914 

9 

|  16.99071 

|  16.54831  | 

|  16.01596 

8 

|  17.41588 

|  16.61207 

|  16.02609 

TabLe  2 
(Sample  Mean) 
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ESPRIT:  Linear  Overlapping  Case 


6  obtained 
from  Bound 


Angle*24° 


9  obtained  from 
chordal  metric 


6  obtained 
from  X 


10  1 

24.26761  | 

24.41945 

|  24.13898 

9  1 

24.41869  | 

24.47776 

|  24.15929 

8  1 

24.68307  | 

24.53942 

|  24.18203 

Angle* 16° 


SNR 

(dB) 

6  obtained 
from  Bound 

6  obtained  from 
chordal  metric 

6  obtgined 
from  X 

10  | 

16.33276 

16.36528 

15.97328 

9  1 

16.52749 

16.41086 

15.97049 

8  1 

16.81625 

16.45882 

15.96841 

TabLe  4 
(Sample  Mean) 
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ESPRIT:  Linear  Non  Overlapping  Case 
&  ESPRIT:  General  Case 


Angle-24® 


SNR 

(dB) 

9  obtained 
from  Bound 

0  obtained  from 
chordal  metric 

9  obtgined 
from  X 

10 

24.51770 

25.58454 

24.04759 

9 

24.82163 

25.94506 

24.03446 

8 

24.30935 

26.38912 

24.01111 

Table  5 
(Sample  Mean) 

Angle-16® 

SNR 

(dB) 

0  obtained 
from  Bound 

9  obtained  from 
chordal  metric 

9  obtained 
from  X 

10 

16.67917 

17.54183 

16.04165 

9 

17.07774 

17.88994 

16.11396 

8 

17.70301 

18.31548 

15.22606 

TabLe  6 
(Sample  Mean) 
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5.3  PERTURBATION  DUB  TO  SENSOR  SPACING 

In  this  case  we  assume  the  environeaent  to  be  noise  free.  However, 
each  sensor  is  assumed  to  be  perturbed  from  its  ideal  position. 

5.3.1  MOVING  WINDOW 

Consider  a  linear  array  of  m  identical  sensors  spaced  a  distance 
D+AD}  where  A Dj  is  the  uncertainty  in  the  spacing  between  the  i-th  and  the 
(i+l)-th  sensors.  Assume  there  are  d  (d<m)  narrowband  sources  located  at 
azimuthal  angles  9^;  k-1,2,.  .  .  d,  which  are  impinging  on  the  array  as 
plane  waves  and  whose  signal  complex  envelopes  are  denoted  by  s^t).  The 
signal  received  at  the  i-th  sensor  is  modeled  as 

d 

^(t.e)-  r  s^DaiCek)  ;  i-1,2 . m  (5. 3. 1-1) 

k-1 

where  denotes  the  response  of  the  perturbed  array  and  is  the 

perturbed  relative  response  of  the  i-th  sensor  to  the  k-th  source.  Note 
that 

ai(e)-a(e)ej^w/c^^i_1^D+aDi^sin^0^ 

-  a(e)ej(i-1)D<“/c>  sin(0)  gj  (tt/c)ADiSin(0)  (5. 3. 1-2) 

where  a(0)  is  the  gain  of  the  sensor  in  the  angular  direction  9.  To  a  first 
order  approximation 

#j(«/c)ADisin(0)  „  i+j(w/c)ADisin(9) 

a  l+j(2itADi/S)sin(0)  (5. 3. 1-3) 

where  $  is  the  wavelength  of  the  signal  wavefront.  Thus,  aj(9)  can  be 
written  as 

ai(9)  -  a(0)e^i-1>D(w/c>  sin(9) 

♦  j ( 2itADi/S)ej (i-l)D(u/c)  sin(0)sin( e)a( 6)  t  (5. 3. 1-4) 

For  simplicity,  let  a|(-a(9|c).  Then 


149 


d 

yi(t,9).  E  •ksk(t)(«J<i-1>D<“/c>  sin<®k> 
k-1 

d 

+j(2nADi/&)  E  aksk(t)(eJ/i-1^D(<^/c^sin<ek>sin(9k);  i«l,...,m.  (5. 3. 1-5) 

k-1 

Notice  that  the  first  part  of  equation  (5. 3. 1-5)  is  just  the  unperturbed 
quantity  y^.  Dropping  the  argument  (t,9)  in  equation  (5. 3. 1-5),  it  can  be 
written  as 

yj.  y*  ♦  Ay*  «  y*  +  ej  .  (5. 3. 1-6) 

(m-L+1)  vectors  Yj^  of  length  L  are  then  formed  where  Yj,  is  given  by 

Xn  *  t  Yn  ^n+l  *  *  *  Yn+L-1 
Yjj  can  be  written  as 

In  “!»♦&»■  (5.3. 1-7) 

From  chapter  2  it  is  shown  that  Yj,  can  be  expressed  in  the  form 

Yn-AB*<n-1>S  (5.3. 1-8) 

where  A,  B,  ♦  and  S  are 

A  *  (f*l  3:2  '  •  •  5d  J 

aj  -  (  1  e^i  .  .  .  ej(L-l)+i] 

B  -  diag  {  aj  a2  .  .  •  ad} 

*  -  diag  [  ej*l  eJ*2  .  .  .  eJ+d  J 

♦k-(»A/c)sin(ek)  ,  k-1,  2,  .  .  d  (5. 3. 1-9) 

S  -  {  sj  s2  •  •  •  sd  JT. 

Similarly  ^  can  be  expressed  as 

En  -  j  (2n/6)  (ADJn  ABGfC"-1^,  (5.3.1-10) 

where  A,  B,  ♦  and  S  are  given  by  equation  (5. 3. 1-9)  and  G  and  [&D]n  are 
G  -  diag  {  sin(9j)  sin(02)  .  •  •  sin(9d)} 

(AD]n  -  diag  {  ADn  AD(n+1)  .  .  .  AD(n+L-l)  )• 

The  two  matrices  Mj  and  Nj  are  then  formed  in  the  usual  way  where 
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"l  ■  I  II  12  •  •  •  II  1 

M,  .  I  Tj  Tj  .  .  .  T(ul)  ]. 

Consider  the  matrices  and  E{N^HM^].  An  element  m^  ^  of  E[M^0M^]  is 

«k,h  -  EiikB-ihi 

where  H  denotes  conjugate  transpose  and  E[.]  denotes  expectation.  Utilizing 
equation  (3. 3. 1-7),  m^*,  becomes 


»k,h  «  E{YkH-Yh]-.E[YkH.EhJ-E£EltH.yhJ^E[EltH.Eh] 


Ellk^'Ihl  can  ke  obtained  in  closed  form  as 


(5.3.1-11) 


d  d 


EfIkH-Ihl-  S  E  SpqapqPpqe-^^1)^ 
q-1  p-1 


where 


^pq-EfsJsp] 

apq*aqap 


(5.3.1-12) 

(5.3.1-13) 

(5.3.1-14) 


Fpq  -  I  eJ(i-D(V  *Q>  . 

H  i-1 


(5.3.1-15) 


Similarly,  EIT^EEj,)!  can  be  expressed  as 


d  d 


ElIkH-Eh]-j(2a/&)  I  I  {  Sp^F^e-W-l)  *q 
q-1  p-1 


x  «J(h-l)*p  sin(9_)} 


where 


h+L-1 

Fpq  -  E  ADj  ej(i*h)(*p"  *q> 


i-h 


(5.3.1-16) 


(5.3.1-17) 


In  a  similar  manner,  EfEj^Yj,)!  is 
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(5.3.1-18) 


d  d 

EIEkH.Yh]-j(2n/S>  £  £  {  SpqapqFjq«-3<k-1)*q 

q-i  p-i 

x  e3(k-l)4p  sin(0q)} 


where 


PQ 


k+L-1 

m  £  ad |  ^q)  • 

i-k 


(5.3.1-19) 


Finally,  following  the  sane  approach,  E| j^°Ep) ]  can  be  written  as 


d  d 

E[EkH.Eh]-j(2lt/S)2  £  £{SpqapqFpq  e-J<k“1>*q 

q-i  p-i 

x  e3(**-l)+p  sin(9q)sin(0p)}  (5.3.1-20) 

where  1 

L-l 

Fj&  -  £  dDr+h  ADr+k  ej(r-k><V  *q>  .  (5.3.1-21) 

r«0 

Note  that  the  matrices  E[M]®M]J  and  E[Ni®M^]  can  be  written  as 
E{MiHMi]  -  M  ♦  E 
and 

EtNj^j)  .  N  ♦  F, 

where  the  kh-th  element  of  N  and  N  are  given  by  EI!kH*IhJ  and  —  L Xk.^lH - IhJ  > 
respectively.  The  kh-th  element  of  E  and  F  are  given  by 

EtIkH-lhI+EI§kH-Ihl+EfEkH-lhl  “d  Ellk^^-Ihl^flk^^-IhJ^Hk.^-IhJ. 
respectively.  Recall  that  ||E||  and  ||F||  are  given  by  the  largest  eigen¬ 
values  of  EHE  and  FBF,  respectively.  Note  that  this  eigenvalue  is  always 
less  than  or  equal  to  the  sum  of  all  the  eigenvalues.  It  is  known  that  the 
the  sun  of  all  eigenvalues  of  a  matrix  is  equal  to  the  trace  of  this 
matrix.  Recall  that  the  matrix  E  is  (m-l)x(m-l).  Let  ej^j,  the  ^h_th  ele- 
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sent  of  E.  It  can  be  shown  that 

(m-L)  (m-L) 

tr(EHE)  -  E  E  hi.il2. 
k-1  h-1 

Note  that  the  square  root  of  tr(D)  is  just  the  definition  of  the  Frobenius 
norm  of  E  defined  as 


f  (m-L)  (m-L) 

|E||f»  Z  1  iei  *  j 

l  k-1  h-1 


(5.3.1-24) 


The  Ilk1*1  element  of  E  is  of  the  the  form 


d  d 


eh.k-K2*/*)  1  ^Spqapqe-J^-D+q  e^-D+p} 
p-1  q-1 

x{F|»qSin(ep)+F£qsin(0q)  +  (j2Ji/S)Fh{;  sin(9p)  sin(9p)}  (5.3.1-25) 


Recall  that  if  a-bc,  then  |a|<|b|+|cj.  Thus, 


d  d 

*  |Spqll>MUI^M^I*(2n/5)|Fhk|)  (5.3.1-26) 
p-1  q-1 


It  is  easy  to  see  that 


ph 

Pd 


L  |  ADBax  |  ;  h#l 

(L-l)|  |  ;  h-1 


(5.3.3-27) 


where  dDj  is  assumed  to  be  zero.  Similarly,  we  can  show  that 


F*  < 
pq- 


L  |  AD, 


max 


( L—  1 )  I  I 


;  k#l 
;  k-1 


(5.3.1-28) 


and 
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(5.3.1-29) 


'  L  |  ADmax  |2  ;  h*l  and  k#l 

<  <  (L-i)  |  ADmax  |2  ;  h-1  or  k«l 
.  (L-2)  |  ADmax  |2  ;  h-1  and  k-1. 

Because  ( 1 ADmay |/S)2  is  very  small,  is  negligible.  Let  R  be  the 
quantity 

d  d 

R  -<2n/S)  I  I  jSpq||apq|  .  (5.3.1-30) 

p-i  q-i 

Ue  can  then  show  that 

'  2L  |  dDmax  | 2  ;  h*l  and  k*l 

lek,hl  <R  1  <2L-D  l^maxl2  t  h"1  or  k«J  (5.3.1-31) 

.  2(L-1)  | ADjnax |2  ;  h-1  and  k-1. 

Using  equation  (5.3.1-24),  it  can  be  shown  that 
1  I E  |  | 2  <  R2  |ADmax|2  [4(L-l)2+2(m-L-l)(?L-l)2+(m-L-l)24L2 ] . 

Similarly,  it  can  be  shown  that 

I  I F  1 1 2  <  R2  I^Wl2  [  (m-L)(2L-l)2+(m-L)(m-L+l)24L2J  • 

Therefore, 


<R  lADaaxI  {[4(L-l)2+2(m-L-l)(2L-i)2+(m-L-l)24L2] 

♦[ (m-L)(2L-l )2+(m-L)(m-L+l )24L2] 

After  some  simplifications,  it  can  be  shown  that 

C<R  lADnaxI  (8L2(m-L)2+(m-L)  (3-12D  +  2)1*.  (5.3.1-34) 

Let  Kl-  [8L2(a-L)2+(m-L) (3-12L)+2 ]^.  The  bound  on  the  chordal  metric  then 
becomes 


(5.3.1-32) 

(5.3.1-33) 
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x(Xi,xi)  i 


(5,3,1-35) 


d  d 

2n(|ADmax|/S)Kl  E  E  |Spq||apq| 
p- 1  q« 1 

\  ♦  (jr^Nx^)2 

Note  that  the  chordal  metric  is  proportional  to  the  correlation  between  the 
sources . 

5.3.2  ESPRIT 

5. 3. 2.1  General  Case 

Consider  a  planar  array  which  consists  of  m  matched  sensor 
doublets  vhose  elements  are  translationally  separated  by  a  displacement 
D+ADj  (  ADj  is  the  uncertainty  at  one  of  the  sensors  of  the  i-th  doublet). 
The  signals  received  at  the  i-th  doublet  are 

d 

*i(t,e)«  E  sk(t)ai(0k) 
k-1 

(5.3.2. 1-1) 
d 

y^t,©)-  E  s».tt)ai(\)  («/c)((  i-l)D*  AOj  )sin(6)  )  _ 

k-1 

To  a  first  order  approximation 

e.1  («/  c)  AD^sin(6)  B  1+.  j  ( y/c )  s  in(  ©) 

•  l*j(2nAD1/5)sin(6),  (5.3.2. 1-2) 

where  S  is  the  vavelength  of  the  signal  wavefront.  Then 

(«/c)( (i-1 )D*ADj )sin( 6) ( i  - 1  )D(<a/c)  sin( 6) 

t 5. 3.2.1  31 

*  sin(6)sin( ©ta(©) . 


It  follows  that 


yi(t,8)-  £  »i<®k>sk<tX«J(i~1)D(W/c)  sin(\) 

k-1 


(5.3,2, 1-4) 


♦  jUittfVS)  £  ai(ek)sk(t)(eJ<i-1>D^w/c)sin<9k)  sin^)  ;  1-1,2, -  .  .,i 
k-1 


Equation  (5.3.2. 1-1)  can  be  written  as 
*i-  *i 

yi-  y{  ♦  Ayj  -  yt  ej  . 


Let  X  and  Y  be  the  vectors 


X  «  {  xj  x2  .  .  .  xffi  }T, 

X  -  {  yi  y2  •  •  •  y®  )T- 

It  is  easy  to  see  that  X  and  Y  can  be  rewritten  as 


(5.3.2. 1-5) 
(5.3.2. 1-6) 


X  -  X 


Y  -  Y  ♦  AY 


(5. 3. 2. 1-7) 
(5.3.2. 1-8) 


Previously,  it  vas  shown  that  X  and  Y  have  the  following  decompositions 


X-Ai§ 

Y-Ai*S. 


(5.3.2. 1-9) 
(5.3.2.1-10) 


yA( t , 8)-  £  a1(6k)sk(t)(e-)(i-1)D^w/c)  *in<®k> 
k-1 


(5. 3. 2. 1-4) 


♦  J(2b6D1/S)  I  ai(ek)sk(t)(eJ<i-1^D(w/c>sin<ek>  sin^)  ;  i-1,2 . m. 

k-1 


Equation  (5.3.2. 1-1)  can  be  vritten  as 


Xj-  Xi 

yi-  yj  *  Ayi  -  Yi  *  *i 


(5.3.2. 1-5) 
(5.3.2. 1-6) 


Let  X  and  Y  be  the  vectors 


(5.3.:. 1-7) 


X  -  {  Xj  X2  .  .  .  Xm  }T, 

!  -  c  yi  h  •  •  •  y*  )T- 

It  is  easy  to  sec  that  X  and  Y  can  be  rewritten  as 

X  -  X 

Y  -  Y  ♦  AY  (5.3.2. 1-8) 

Previously,  it  was  shown  that  X  and  Y  have  the  following  decompositions 

X.A1S  (5. 3. 2. 1-9) 

Y*AX ♦§  (5,3.2.1-10) 

where  A,  ♦  and  S  are 

Aj  -  [a(®i)  a(®2^  •  •  •  a(0<})] 
a(0j)  -  [  aj(0j)  a2(0j)  -  .  .  am(0i>] 

♦  -  diag  [  eJ*l  .  .  .  e^d  ] 

S  •  {  sj  s2  .  .  .  s^  )T. 

Similarly  AY  can  be  expressed  as 

AY  -  j  (2n/5)  ( AD]  AjGiS,  (5.3.2.1-11) 

where  A,  ♦  and  S  have  been  defined  earlier  and  G  and  [AD]  are 
G  «  diag  {  siniOj)  sin(@2)  .  .  .  sin(8d)} 

[AD]  -  diag  {  ADi  AD2  .  .  .  ADm  } • 

Let  M  and  N  be  the  matrices 

M-E[X  XH  ]-E[X  X]-M  (5.3.2.1-12) 

M-E[X  YH  ]-E[X(YH+AYH)1-E[X  YH]+E[X  AYh]-N*AN.  (5.3.2.1-13) 

The  error  matrices  E  and  F  are  given  by 

E-0  (5.3.2.1-14) 

F-AN-E[X  AYH1— j(2it/&)A1S*HGHA1H[AD]H.  (5.3.2.1-15) 

Recall  that 


(5.3.2.1-16) 


is? 


It  is  possible  to  show  that  the  hk-th  element  of  the  matrix  F  is 


*h,k«-j<2K/*>ADhTh,k 


(5.3.2.1-17) 


where 


d  d 


and 


Therefore, 


Th,k  *  ^  ^  ah^®p^ak  (®q)Spq  sin(9q)  e  J^q 

P-1  q-i 


Spq  “  ElSpSq  ]. 


(5.3.2.1-18) 


d  d 


|fhfk|  $(2*/S)  |ADmax|  E  E  |eh(©p)||ak*(eq)||Spq|.  (5.3.2.1-19) 

p-i  q-i 


For  omni-directional  sensors,  a(6)-l,  we  get 


d  d 

I'h.kl  s  1“  max  I  r  E  I spq I ' 

p-l  q-i 


(5.3.2.1-20) 


The  Frobenius  norm  of  F  is  given  by 


l|F||  -  {  ”  "|fh,kl2)‘ 

l  h-1  k-1  J 


(5.3.2.1-21) 


Thus, 


||F||  <  (2ll/S)  |bDmax| 


m  m  d  d  V* 

E  E  E  E  |S  |2  . 

h-1  k-1  p-l  q-1  ) 


which  reduces  to 


(  d  d 

||F||  <  (2a/i)  | ADmax (  m  l  E  (Spp  |2  . 

{  p-l  q-1  J 


(5.3.2.1-22) 


(5.3.2.1-23) 
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The  bound  on  the  chordal  metric  then  becomes 


d  d 

2n(|dDmax|/$)K2  E  I  |Spq| 
p-1  q«l 

X(Xi,Xi)  <  - .  (5.3.2.1-24) 

J  (2iHMXi)2  ♦  (^^Nxj)2 


where 

K2«m. 

5. 3. 2. 2  Linear  Array:  Overlapping  Doublets  Case 

In  this  section  ve  consider  ESFRIT  where  a  linear  array  composed 
of  m  sensors  is  used  to  solve  for  the  angles  of  arrival.  As  for  the  case  of 
the  moving  window,  we  assume  that  the  i-th  sensor  is  displaced  by  an  amount 
ADj  with  respect  to  the  reference  sensor  which  we  assume  as  the  first 
sensor.  As  seen  in  section  5.3.1,  the  received  signal  at  the  i-th  sensor  is 
given  by 

d 

yj(t,8)»  E  a^s^f t)(e^ (i**l)D(«/e)  sinC©^) 
k-1 

♦  (5. 3. 2. 2-1) 

d 

j(2nAD1/S)  E  *ksk(t)(*^*~1^w/C^sin^0^  sinOj,)  ;  i-1,2,.  .  .,m. 
k-1 

Two  arrays  Yj  and  Y2  are  then  formed  from  this  data  where 

11  -  {  yi  Y2  •  •  •  >’m-l  )T 

12  -  (  V2  Y3  *  •  •  ym  )T* 

Recalling  that  y^  can  be  expressed  as 

Yi  “  Yi  +  *i  5  i-1*  2,  •  •  • 

where  yj  is  the  unperturbed  data  and  ej  is  the  error  in  yj.  Therefore  Yj 
and  Yj  are  given  by 
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|i  -  Ti  ♦  Ex. 

12  *  12  ♦  12* 

II  an<*  12  k*v®  the  following  decompositions 

y1.a2bs 

Y  t  »A  2  B  $S 

where  A,  B,  ♦  and  S  are 

A2  -  l£i  a2  •  .  .  ad  ] 

at  -  [  1  *i*i  .  .  .  eJ<m-2>*i] 

B  -  diag  {a!  a2  .  .  .  ad} 

*  -  diag  [  «3*1  .  .  .  el*d  1 


(5. 3. 2. 2-2) 


(5. 3. 2. 2-3) 


(5. 3. 2. 2-4) 


(5. 3. 2. 2-5) 


♦jt»(«A/c)sin(6|t)  ,  Jt-1,  2,  .  .  d.  (5. 3. 2. 2-6) 

Similarly  Ej  and  ^  can  be  expressed  as 

Ex  -  j  (2a/ S)  [ AD]X  A2BGS,  (5. 3. 2. 2-7) 

*2  -  j  (2a/ S)  [ADJ2  A2BG#S,  (5. 3. 2. 2-8) 

where  A,  B,  #  and  S  have  been  defined  earlier  and  G,  [ADJj  and  [AD}2are 
G  -  diag  (  sin(8^)  sini^)  .  .  .  sin(8d)}, 

(ADJi  -  diag  {  ADj  AD2  .  .  .  ADm_1  ), 

[ AD  1 2  -  diag  {  AD2  AD3  .  .  .  ADa  ). 

Two  matrices  M  and  N  are  then  formed  where 
M-E[Yi  Y^I-EKYi  *  IlKIl  ♦  li)Hl, 

N-EIXl  !H2l-El(Il  *  E1)(Y2  «■  E2)H]. 


These  can  be  decomposed  into 

M-EtY!  XB1 1 ^Elll  EHi ]^E[EX  XHi 1^E[Ei  E^j J *M+E,  (5. 3. 2. 2-9) 

N-EIXi  XH2)+EtXl  IH2 1 ^ElEl  XH2l*EtIl  EE2)*N+F,  (5.3.2.2-10) 


where 


H-EIXi  XHll. 
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N-E(Xi  Ihi]+eHi  Ihi1+e(Ii  ehj1, 

E-EIXi  XH2l * 

F-E(Xi  EH2l+EHl  IH2l+EJIl  lH2l  * 

To  obtain  a  bound  on  the  chordal  metric  one  needs  to  explicitly  express  c 
in  terms  of  the  parameters  of  interest.  One  has  thus  to  get  the  Euclidean 
norms  of  E  and  F.  It  can  be  shown  that 

E l II  IH1 1 — j ( 2 R/ 5) A2BSBHGHA2H t ADI iH, 

E[Ei  XHlI*EtIi 

E[§i  EH1]-(2a/5)2  [ AD]1A2GBSBEGEA2E[ ADJ^®, 

E(Xl  EH2I— j(2K/5)A2BSBHGH*A2H[AD]1H, 

E(Ei  XH2 1-J (2«/S) [ AD] 1A2GBSBh*A2h, 

Et|i  EH2]-(2ll/S)2  [AD]1A2GBSBHGH*A2H[AD]2H. 

Consider  the  matrix  E]Yi  EEj].  The  hk-th  element  of  this  matrix  can  be  ex¬ 
pressed  as 

-j(2it/i)ADh  Thk, 

vhere 

d  d 

Thk  -  I  I  apq  Spq  sin(0q)  e-j(k-l)*q, 

p-i  q-i 

and  apq  and  Spq  have  been  defined  earlier.  Since  the  matrix  ElEj  XHl )  vas 
shovn  to  be  E [ X^ll^*  the  hk-th  element  of  this  matrix  can  be  written  as 
j(2lt/S)ADh  T^*. 

Because  the  matrix  E[Ej  E*^]  is  given  by  (2il/6)2  [  AD]  ^AGBSB^G^^A^]  AD]  a 
term  of  the  form  (AD/5)2  will  be  contained  in  all  the  elements  of  this 
matrix  causing  the  matrix  to  be  negligible.  Therefore,  the  hk-th  element  of 
the  matrix  E  is  gi  ven  by 

«hk  "  -j(2»/8)ADh  Thk  *  j(2tl/5)ADh  T^*,  (5.3.2.2-11) 


lbl 


(5.3.2.2-12) 


ehk  «  j(2«/«)ADh  (“Thk  *  Tkh*>* 

Therefore, 

d  d 

|ehk |  <2(2n/S)  |ADmax|  E  E  |apq  |  |Spq|  (5.3.2.2-13) 

P-1  q-1 

The  Frobenius  norm  of  E  is  given  by 

(mm 

l|E||  -  E  E  |ehk|2  (5.3.2.2-14) 

l  h-1  k-1  ) 


Thus 


d  d 

||E||  <  <2ll/A)  |ADmax|[2(m-2)(m-l)]*  E  £  l^pqllS^I  .  (5.3.2.2-15) 

P-1  q-1 

The  second  step  is  to  compute  ||F||.  For  this  recall  that 
F-E[Yi  Eh21+E[E1  Yh2}^E[E1  EH2], 

where 

Edl  EH21— j(2lt/S)A2BSBHGH*A2HlADl1H, 

E(Ej  Y®2  ]-j  (2»l/ 4) [  AD]  iA2GBSBEfA2E, 

E[Ei  EH2]-(2it/&)2  (AD]1A2GBSBHGH*A2H[AD]2H. 

Because  E[Ej  EH2]  will  have  a  factor  of  (AD/4)2  which  tends  to  zero,  this 
term  becomes  negligible  and  is  omitted  in  the  computation  of  the  matrix  F. 
An  element  hk-th  of  the  matrix  E[Y^  EH2]  is  of  the  form 

d  d 

-j (2 a/ 4) ADj,  E  E  apq  Spq  sin(6q)  e^h-D+p  e~ik^q. 

P-1  q-1 

The  hk-th  element  of  the  matrix  E[Ej  Y®2)  is  of  the  form 
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d  d 

r  E  apq  Spq  sin(0p)  e^'1^  e'^^q. 

P-1  q-1 


The  hk-th  element  of  the  matrix  F,  fj^,  is  given  by 


d  d 

fhk-j( 211/5)^  E  E  apqSpq  ei  (h~l)  +?e~)k+q(sin(  ep)+sin(0q) ) .  (5.3.2.2-16) 
p-1  q«l 


Therefore, 


d  d 

|fhk!  <  (2it/&)|fiDaax|  E  E  |&pq  |  |Spq|,  (5.3.2.2-17) 

p-1  q-1 


and 


d  d 

|F|  <  2it(|bDmax|/6)I(m-l)2+(m-l)(m-2)]%  E  E  |apq  |  |Spq  | .  (5.3.2.2-18) 

p-1  q-1 


Recall 


«-J  INI2  •  IlFll2' 

d  d 

e<2n( | ADmax |/&) [(m-l)^+(m-l)(m-2)+2(m-2)(m-l)]^  E  E  |apq  |  |Spq|. 

p-1  q-1 

d  d 

e  <  2n( |dDBax|/6)[(m-l)(4m-7)]lfe  E  E  |apq |  [Spq|.  (5.3.2.2-19) 

p-1  q-1 


Let  R3  -  ( (m-l)(4m-7) ]^.  The  bound  on  the  chordal  metric  then  becomes 


d  d 

2n(|ADmax|/6)K3  E  E  |Spq|!apq| 
p-1  q-1 

X(Xi,Xi)  <  - .  (5.3.2.2-20) 

J  (yiHMxi)2  *  (£iHNxi)2 
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S.3.2.3  Linear  Array:  Non-Overlapping  Doublets  Case 

In  this  section  ve  assume  that  tvo  adjacent  sensors  form  a  pair  so 
that  tvo  non-ovelapping  arrays  Y3  and  Y2  are  formed  where 

11  -  {  yi  y3  •  •  •  yB-i  )T 

12  •  t  Y2  Y4  •  ■  •  ym  )T- 

m  is  assumed  to  be  a  multiple  of  2.  Recalling  that  yj  can  be  expressed  as 
Yi  *  yj  ♦  e^  ;  i“l»  2,  >,m, 

where  y*  is  the  unperturbed  data  and  e^  is  the  error  in  y±.  Therefore,  Y3 
and  Y2  can  be  expressed  as 


II  “  II  -  Iv 

h  - 12  ♦  i2- 

Yj  and  Y2  have  the  following  decomposition 

Yj^BS  (5. 3.2. 3-1) 

Y2-A3B*S  (5. 3. 2. 3-2) 

vhere  A3,  B,  ♦  and  S  are 

A3  *  til  i2  •  •  '  5d  J 

aj  -  [  1  eJ2*i  .  .  .  ej("-2>*i] 

B  «  diag  {  aj  82-..  a^) 

♦  «  diag  [  eJ^l  e^^2  .  .  .  e^d  ] 


and 


tfe-CcoA/Osinfe^)  ,  k«l,  2,  ....  d. 

Similarly  E3  and  Ej  can  be  written  as 

Ej  -  j  (2JI/S)  IADI3  A3BGS,  (5. 3. 2. 3-3) 

E2  -  j  ( 2 it/ 6)  [AD]2  A3BG*S,  (5. 3. 2. 3-4) 

vhere  A3,  B,  ♦  and  S  have  been  defined  earlier  and  G.  [AD] 3  and  [AD)2are 
G  -  diag  {  sin(©3)  sin(82)  •  •  •  sinie^)), 

[AD]3  -  diag  {  AD 3  AD3  .  .  .  ADB_3  }, 
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[AD] 2  -  diag  {  AD2  AD4  .  .  .  ADm  } . 

Tvo  matrices  H  and  N  are  then  formed  where 
N-EIX!  Y^l-EUYi  +  IiHIl  +  El)H], 

N-ElYj  |H2]»E[(Y1  +  E1)(Y2  +  Ej)8]. 

These  can  be  decomposed  into 

M-ElY!  YH1]-t-E[Y1  EH1]fE[E1  YH1]+E[E1  E^l-M+E,  (5. 3. 2. 3-5) 

N-EIY!  YH2]-*-E[Y1  EH2]+E[E1  YH2]+E[E1  EH2]-N+F,  (5. 3. 2. 3-6) 

vhere 

M.E[Yi  IHll» 

H-ElTi  ^il^EIE!  Y^l+EtEj  EH2  ] , 

E-E[Yi  YH2], 

F-E[Yi  E®2]+E[E1  Yh2]+E[E1  Eh2], 

It  can  be  shovn  that 

E III  IH1 1  — J ( 2 */  S) A3BSBHGHA3H [ AD ] j8 , 

Elll  Ihi1*e(ei  IHi]H, 

EHl  EH1]-(2n/S)2  [ADJ1A3GBSBaGHA3H[AD]1H, 

EHl  EH21— j(2it/5)A3BSBHGH*A3H[AD]1H, 

ElEl  YH2]-j(2lt/S)[AD]1A3GBSBH*A3H, 

E [ Ei  EH2]-(2lt/5)2  [AD]1A3GBSBHGH*A3H[AD]2H. 

Consider  the  matrix  E[Yj  E8jJ  .  The  hk-th  element  of  this  matrix  can  be  ex¬ 
pressed  as 

-j(2Jt/6)AD2h_1  Thk 

vhere 

d  d 

Thk  -  I  I  apq  Spq  sin(6q)  e:K2h-2>*p  e"j <2~2>^q;  h , k-1 . 2 , . . . , (m/2) 
P-1  q-1 

and  Spq  and  Spq  have  been  defined  earlier.  Since  the  matrix  E [ E3  vas 
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shown  to  be  E[Ej  Y® i ] ® »  the  hk-th  element  of  this  matrix  can  be  written  as 
j(2n/8)bD2h_1  Thk*. 

Because  there  will  be  a  term  (&D/$>2  in  the  matrix  E[E^  E®j]  which  tends  to 
zero,  the  hk-th  element  of  the  matrix  E  is  approximated  by 

«hk  *  -j(2n/&)bD^  +  j(2Jt/S)bdj1  Tj^*,  (5. 3. 2. 3-7) 

«hk  *  j(2«/4)ADh  (-Thk  ♦  T^*).  (5. 3. 2. 3-8) 

Therefore, 

d  d 

lehk I  <  2(2lt/&)  \tomax\  Z  Z  |ap  |  |Spq|. 

p-i  q-i 

The  Frobenius  norm  of  F  is  given  by 

(  m/2  m/2 

||E||  .II  |ehk|2  . 

I  h.l  k.l  ) 

Thus 

d  d 

||E||  <  (2n/S)  |ADmax|[2(m/2-l)(m/2)]%  E  Z  japq  j  jspq  | .  (5.3.2.3-11) 

p-i  q-i 

The  second  step  is  to  compute  ||F||.  Recall  that 
F-Etll  E^l+EIE!  Yh2]+E[Ei  Eh2]. 

The  element  E[Ej  E®2]  will  be  neglected  because  it  will  have  a  factor  of 
(AD/&)2  which  tends  to  zero.  The  hk-th  element  of  the  matrix  E [ EH2)  is 
of  the  form 

d  d 

-j(2it/$)4D2h  Z  Z  apq  Spq  sin(9q)  ej(2h"2>*p  e-j(2k-l)*q. 

p»l  q*l 

The  hk-th  element  of  the  matrix  E[Ej  Y®2)  is  of  the  form 


(5. 3. 2. 3-9) 


(5.3.2.3-10) 


d  d 

j(2*/S)AD2h_1  E  Z  apq  Spq  sin(6p)  e^2*1'1)*?  e-j (2k~2> *q . 

P-1  q-1 

Therefore, 

d  d 

IfhiJ  <  (2H/S)  |ADmax|  r  Z  |apq|  |spq|,  (5.3.2.3-12) 

p-i  q-i 


d  d 

|F|  <  2n(  |dDoax|/5)t(m/2)2+(m/2-l)(m/2)]!,i  Z  Z  |apq  |  |Spq|.  (5.3.2.3-13) 

P-1  q-1 

Recall 


d  d 

e<2n(  |ADmax|/6M(m/2)2+(m/2-lHm/2)+2(m/2)(m/2-l)]%  E  E  |apq  |  |Spq|. 


d  d 

e  <  2rc( | ADmax |/S)[(m/2)(2m-3)]*  E  E|apq|jSpq|.  (5.3.2.3-14) 

p-1  q«l 


Let  K.4  *  [ (m/2)(2m-3) ]^.  The  bound  on  the  chordal  metric  then  becomes 


XC  ) 


d  d 

2*(|ADmaxl/5>K4  2  2  |Spq||apq| 


(5.3.2.3-15) 


5.3.3  COMPUTER  SIMULATION 

In  this  section,  we  studied  the  effects  of  errors  due  to  sensor 
spacing  on  the  performance  of  the  3  algorithms  discussed  earlier.  The  model 
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used  consisted  of  tvo  incoherent  (d«2)  incident  on  a  linear  array  consist¬ 
ing  of  eight  sensors  (m*8).  The  sources  are  assumed  to  be  located  at  0^=16° 
and  ©2-24®.  For  simplicity,  the  case  of  omnidirectional  sensors  vas  assumed 
for  all  three  cases.  In  the  simulation  the  case  of  perfect  sensor  spacing 
vas  first  considered.  100  snapshots  vere  used  to  obtain  the  matrices  M  and 
N.  The  process  vas  repeated  50  times  and  the  results  averaged  to  obtain 
nominal  values  for  Xj ,  and  ;i«l,2.  A  random  perturbation  vith  a  maxi¬ 
mum  AD  varying  from  D/100  to  D/1000  vas  then  introduced  and  the  procedure 
used  in  the  unperturbed  case  vas  repeated.  D  vas  assumed  to  be  equal  to 
half  the  vavelength  so  that  a>4/c*n  (4  being  the  vavelength).  The  computed 
results  arc  shovn  in  tables  5.7  to  5.14.  If  the  error  is  small  enough,  then 
the  bounds  derived  in  this  section  give  acceptable  results.  Hovever,  if  the 
error  large,  then  the  conditions  for  vhich  these  bounds  have  beer  derived 
do  not  hold  any  more  and  therefore  the  bounds  are  not  applicable. 
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Moving  Vindev 


Angle-24® 


AD 

0  obtained 
from  Bound 

8  obtained  from 
chordal  metric 

6  obtained 
from  X 

0/100 

26,4719a 

24.00648 

24.00001 

D/300 

24.81639 

24.00216 

24,00001 

D/500 

24.48902 

24.00130 

24 . 00000 

D '  700 

24.34906 

24.00093 

24.00000 

0  1000 

24.24422 

24.00064 

24 . 00000 

Table  5.7 
(Sample  Mean) 


Angle- 16® 


AD 

0  obtained 
from  Bound 

8  obtained  from 
chordal  metric 

8  obtained 
from  \ 

D/100 

20.68584 

16.00554 

16. 000 '0 

D  300  | 

17.53180 

16.00184  j 

1  P.0002 4 

D  500  | 

16. <31655 

16.00110 

IP.  OOOU 

D-'  700  | 

16 . 6530'  j 

16.00070 

IP. 00010  1 

D  1000 

16. 45 '43 

16.00055 

1 p . 0000 ' 

Table  5 . 8 
(Sample  Mean' 


>*> 


ESPRIT:  Linear  Overlapping  Case 


Angle-24' 


AD 


6  obtained 
from  Bound 


8  obtained  from 
chordal  metric 


8  obtained 
from  X 


D/100  j  25.67951  |  24.00629 


24.00038 


D/300 


24.55655 


24.00210 


24.00015 


D/500  |  24.33357  |  24.00126 


24.00009 


D/700  |  24.23815 


24.00090 


24.00006 


0/ 1000  I  24.16665  |  24.00063 


24.00005 


Table  5.9 
(Sample  Mean) 


Angle-16® 


AD 

!  6  obtained 
from  Bound 

8  obtained  from 
chordal  metric 

■v 

0  obtained 
from  X 

D/100 

18.17524 

16.00500 

1 6 . 000 ' 7 

D/300 

16.72113 

16.00167 

16.00021 

D/500 

16.43227 

16.00100 

16.00013 

D/700 

16.30865 

16.00071 

16.00009 

D/1000 

16. 215°9 

16.00050 

16.000(76 

TabLe  5.10 
(Sample  Mean! 


ESPRIT:  Linear  Non  Overlapping  Case 


Angle-24® 


AO 

9  obtained 
from  Bound 

9  obtained  from 
chordal  metric 

9  obtained 
from  X 

P/100 

25.76599 

24.01485 

|  24.00153 

D/300 

24.58539 

24.00493 

j  24.00053 

D/500 

24 . 35086 

24.00296 

|  24.00032 

D/700 

24.25051 

24.00212 

|  24.00023 

D/IOOO 

24.17530 

74.00148 

|  24.00016 

Table  5.11 
(Sample  Mean) 


Angle-16® 


AD 

9  obtained 
from  Bound 

9  obtained  from 
chordal  metric 

9  obtained 
from  X 

D/100 

18.25734 

16.01279 

16.00038 

D/300 

|  16.74875 

16.00426 

16.00012 

D/500 

|  16.44886 

16.00255 

16.00007 

D/700 

[  16.32050 

16.00183 

16.00005 

D/1000 

16.00128 

16 . 00003 

TabLe  5.12 
(Sample  Mean) 
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ESrRIT:  General  Case 


Angle»24° 

AO 

6  obtained 

8  obtained  from 

6  obt§ined 

from  Bound 

chordal  metric 

from  X 

D/100 

25.57809 

24.01101 

23.99986 

D/300 

24.52344 

24.00366 

23.99995 

D/500 

24.31377 

24.00220 

23.99997 

D/700 

24.22404 

24.00157 

23.99998 

D/IOOO 

24.15678 

24.00110 

23.99998 

Table  5.13 
(Sample  Mean) 


Angle«16® 


AD 

9  obtained 
from  Bound 

6  obtained  from 
chordal  metric 

9  obtained 
from  X 

D/100 

18-01733 

16.01102 

15 . 99943 

D/300 

16.66954 

16.00368 

15.99981 

D/500 

16.40142 

16.00221 

15.99989 

D/700 

16.28663 

16.00158 

15.^9992 

D/1000 

16.20060 

16.00111 

15.00004 

TabLe  5.14 
(Sample  Mean) 


CHAPTER  6 


CONCLUSION  AND  FUTURE  RESEARCH 


6.1  CONCLUSION 

In  this  chapter  ve  discuss  the  contributions  of  this  work 
which  dealt  exclusively  vith  the  non  search  procedure  known  as  the  Matrix 
Pencil  Approach.  ESPRIT  and  the  Moving  Window  are  but  two  of  the  operators 
that  can  be  used  in  the  formulation  of  this  approach  which  is  based  on  the 
generalized  eigenvalue  decomposition  of  two  matrices  generated  from  the 
received  data.  Special  attention  is  given  to  the  Moving  Window  since  it  was 
shown  to  apply  even  in  the  case  of  fully  correlate'*  signals.  However,  the 
method  as  applied  in  [11  did  not  perform  as  was  expected  especially  in 
cases  of  low  signal  to  noise  ratio.  We  have  shown  that  the  separation  of 
the  signal  and  noise  subspaces  is  possible  with  the  use  of  a  window  of 
length  L  with  d<L<(m-d)  where  m  is  the  number  of  sensors  and  d  is  the  num¬ 
ber  of  sources.  This,  in  turn,  allowed  us  to  consider  only  those  eigen¬ 
values  which  are  related  to  the  signal  subspace  by  using  a  singular  value 

decomposition  (SVD).  In  [11  a  window  of  length  L-m-d  was  recommended.  We 
have  explained  why  this  choice  is  actually  the  worst  one  since  it  results 
in  matrices  of  dimension  (dxd)  which  do  not  permit  recognition  of  the  tvo 
subspaces. 

Most  of  the  proposed  high  resolution  techniques  in  direction  find¬ 
ing  treat  each  sensor  in  the  array  as  if  it  exists  by  itself.  In  practice 

however,  mutual  coupling  exists  and  is  very  strong  if  the  separation  be¬ 

tween  the  sensor  is  small.  This  can  significantly  alter  the  structure  of 
the  matrices  involved  in  the  formulation  of  the  proposed  algorithms  vhich. 
in  turn,  drastically  degrades  their  performance.  In  chapter  3  ve  have  pro- 
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posed  a  model  vhich  takes  into  account  the  effects  of  mutual  coupling  be¬ 
tween  the  sensor  elements  of  the  array.  Under  these  conditions,  ve  have 
studied  the  performance  of  the  Matrix  Pencil  Approach  and  ve  have  shown 
that  the  decompositions  needed  in  this  formulation  are  not  possible.  Ve 
proposed  tvo  methods  to  solve  this  problem.  The  first  method  consists  of 
obtaining  an  estimate  of  the  incident  signal  vector.  A  minimum  mean-squared 
error  estimation  was  then  performed  and  an  estimate  was  found.  It  was 
shovn,  through  computer  simulations,  that  the  angles  of  arrival  of  the 
sources  are  well  estimated  using  this  estimate.  In  the  second  method,  some 
pre-processing  was  needed  in  order  to  generate  the  desired  incident  signal 
vectors.  This  vas  referred  to  as  the  Direct  Method.  Several  schemes  have 
been  proposed  depending  on  the  nature  of  the  algorithm  used.  All  schemes 
have  been  shovn  to  be  successful  in  estimating  the  angular  locations  of  the 
sources. 

The  previous  analysis  dealt  with  narrowband  signals.  The  modeling 
used  there  is  not  appropriate  when  dealing  with  wideband  sources.  Chapter  * 
deals  with  signals  of  this  nature.  Ve  have  devised  three  techniques  the 
first  of  vhich  is  original  in  the  sense  that  the  Matrix  Pencil  Approach  is 
utilized  with  a  signal  model  not  used  previously  in  other  approaches.  The 
signals  are  identified  by  their  natural  frequencies  and  their  angles  of  ar¬ 
rival.  This  modeling  is  appropriate  vhen  the  source  signals  are  non  sta¬ 
tionary.  Three  matrix  pencils  have  been  generated  from  the  data.  The  rank 
reducing  values  of  the  first  matrix  pencil  allows  us  to  generate  estimates 
of  the  natural  frequencies.  The  rank  reducing  values  of  a  second  matrix 
pencil  are  shovn  to  be  related  to  both  the  angles  of  arrival  and  the  natu¬ 
ral  frequencies  of  the  sources.  At  this  stage,  it  is  not  apparent  vhich 
natural  frequencies  go  with  vhich  angles  of  arrival.  The  rank  reducing 
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values  of  a  third  matrix  pencil  are  used  to  eliminate  any  ambiguities  that 
could  arise.  The  second  method  utilizes  the  same  model  used  by  Su  and  Morf. 
However,  the  array  configuration  used  is  similar  to  the  first  method.  The 
sources  are  assumed  to  be  linear  systems  driven  by  white  noise  sequences.  A 
scheme  was  devised  in  which  the  angles  of  arrival  could  be  solved  for  with 
the  knowledge  of  the  system  poles.  These  poles  are  shown  to  be  a  mixture  of 
the  source  poles  and  the  sensor  poles.  The  analysis  is  carried  out  on  the 
unit  circle  by  using  a  discrete  Fourier  transform  on  the  data  sequences. 

The  third  method  makes  use  of  the  CSS  of  Vang  and  Kaveh.  This  method  was 
used  in  conjunction  with  ESPRIT  and  the  Moving  Window.  The  method  is  shown 
to  perform  very  well  and  the  sample  variances  of  the  angle  estimates  are 
shown  to  closely  follow  the  Cramer-Rao  Lower  Bound  (CRLB). 

Chapter  5  analyzes  the  effects  of  the  noise  and  perturbations  due 
to  sensor  spacing  on  the  performance  of  ESPRIT  and  the  Moving  Window.  A 
measure,  termed  the  chordal  metric,  was  introduced.  The  chordal  metric  is 
shown  to  be  a  function  of  the  true  and  perturbed  angles  of  arrival.  Geom¬ 
etric  upper  bounds  have  been  derived  for  the  Moving  Window  and  ESPRIT  oper¬ 
ators.  The  proposed  bounds  give  insight  into  performance  degradation  when 
ideal  modeling  is  not  met. 

6.2  PUTURI  WORK 

The  Matrix  Pencil  Approach  is  based  on  exact  knowledge  of  the  num¬ 
ber  of  sources.  Several  methods  have  been  proposed  in  the  case  of  Gaussian 
signals  [40,41].  These  methods  are  shown  to  be  very  effective  with  respect 
to  some  set  criteria.  Special  efforts  should  be  devoted  to  non  Gaussian 
cases.  Also,  the  wideband  methods  described  in  chapter  *  assume  the  number 
of  natural  frequencies  (first  method)  and  the  number  of  poles  (second  meth- 


od)  are  known.  Significant  distortion  will  arise  if  this  number  is  un¬ 
derestimated.  More  effort  is  needed  for  this  particular  case. 


Ve  have  studied  the  effects  of  mutual  coupling  between  the  array 
elements  in  the  narrowband  case  where  only  a  single  carrier  frequency  is 
assumed.  In  the  case  of  wideband  sources,  the  mutual  impedances  become  fre¬ 
quency  dependent.  This  significantly  changes  the  nature  of  the  signal 
modeling.  The  method  that  we  have  used  should  be  generalized  to  the  case  of 
wideband  signals. 

Notice  also  that  the  bounds  derived  for  the  chordal  metric  in 
chapter  5  are  not  very  tight.  This  is  mainly  due  to  the  procedures  used  in 
evaluating  the  Frobenius  norms  of  the  error  matrices.  There  exist  other 
techniques  to  evaluate  these  norms.  One  can  certainly  tighten  these  bounds 
in  order  to  obtain  more  insight  into  performance  degradation  when  ideal 
conditions  are  not  met. 

An  interesting  case  arises  when  one  mounts  an  array  of  sensors  on 
an  airplane.  The  vibrations  of  the  airplane  will  cause  the  sensors  to  be 
displaced  from  their  ideal  positions.  A  two  dimensional  perturbation  analy¬ 
sis  is  needed  to  evaluate  the  chordal  metric.  Also,  one  is  expected  to 
study  the  effects  of  the  structure  on  the  mutual  impedances. 

Finally  B.  Ouibrahim  (1J  proposed  a  third  operator,  called  the 
summation  operator,  that  can  be  used  in  the  formulation  of  the  matrix  pen¬ 
cil.  The  work  that  we  have  developed  here  can  be  easily  generalized  using 
this  operator. 


APPENDIX 


COMPUTATION  OF  THE  CRLB 

Consider  a  linear  uniformly  spaced  array  consisting  of  m  wideband 
sensors  and  let  there  be  d  wideband  sources  (d<m)  located  in  the  far  field 
and  emitting  signals  arriving  at  the  array  from  direction  6^;  i«l,2,...,d. 
The  observed  data  vector  X  at  a  frequency  can  be  expressed  as 

XC®!)  «  A^Sf^)  ♦  N^)  ;  1-1,  2,  .  .  .  ,  L  (A.l) 

where 

A  is  the  direction  matrix 
S  is  the  source  vector 

and 

N  is  the  additive  noise  vector. 

Assume  the  noise  components  to  be  statistically  independent  zero 
mean  random  variables  with  variance  a2 -  Assume  also  that  S  is  a  zero  mean 
random  vector.  Let  R  be  the  covariance  matrix  of  the  observed  data  vector 
X.  Let  8  denote  a  parameter  vector  whose  elements  consist  of  the  angles  of 
arrival  and  statistical  parameters  related  to  the  signal  and  noise  complex 
envelopes.  The  joint  probability  density  function  of  X  given  8  is  given  by 
f(X/8).  (2it)-(®/2>  {det(R)}_lfe  exp{-(l/2)XHR_1X  }.  (A. 2) 

Therefore, 

Log{f(X/8)}-  -(m/2)Log(2n)-(l/2)Log(det(R)>  -( 1  2 >XHR~1X.  (A. 3) 
Taking  into  account  all  the  frequency  components  and  assuming  statistical 
independence  from  one  band  to  the  next,  we  obtain 

L  L 

Lcg(f(X/8)}-  C  -(1/2)  r  Logtdet (R) )  -(1  2)  Z  XfiR_1X.  tA.-) 

1-1  1-1 
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where  C  is  a  constant.  Let  9  be  an  unbiased  estimator  of  8.  It  is  known 
that 

Var(8)  >  J-1(0>  (A. 5) 

where  J  is  the  Fisher  information  matrix  whose  ij-th  entry  is 

"  E  t  (  3Log(  f  (  X/8)  )  /  30J  )  (  3Log(  f  ( X/ 8)  ) /  30j  ) } .  (A. 6) 

It  has  been  shown  (89]  that 

I J  ( ©)  J  i ,  j  -  (1/2)  Tr  {  (R-1  aR/SejXR-1  3R/38.j)}  (A. 7) 

where  Tr(B)  denotes  the  trace  of  the  matrix  B. 

For  the  sake  of  clarity,  assume  that  2  correlated  sources  s^  and 
S2  impinge  on  the  array  from  directions  8^  and  respectively.  Let  p  be 
their  correlation  coefficient.  Assuming  that  the  noise  components  are  inde¬ 
pendent  2ero  mean  random  variables  with  variance  on2 ,  the  covariance  matrix 
R  can  be  written  as 

R  -  E  [X  XH]  -  A  S  AS  «■  4  1  (A.S ) 

where  I  is  the  (mxm)  identity  matrix.  The  matrix  S  can  be  expressed  as 

2 

o  p  0^02 
p  02  °2  ^ 

where  the  variances  of  s^  and  S2  are  denoted  by  oj^  and  o2S  respectively. 
The  signal  to  noise  ratios  SNR}  and  SNR2  are  defined  as 
SNR^  -  101og{  oj2/^}, 

SNR2  -  101og{  o22/o2). 

Therefore,  oj2  and  o22  are  given  by 
O}2  .  o2  iotSNRi/10), 
o22  -  o2  io(SN^2/10>. 


The  matrix  S  can  be  rewritten  as 


'  s2  i0(SNR!/10)  p  gl  10(SNR1+SNR2)/20) 

s  » 

p  ^  10((SNR1+SNR2)/20)  g2  10(SNR2/10) 
The  matrix  A  is  of  the  form 


A 


'  1  1 

eJ  *1  *2 

•  • 

*  r 

ej(m-l)4i  ej(m-l)42  J  f 


where  *  (®A/c)sin(6x ) .  The  parameter  vector  8  is  given  by 
&  «  [  ©j,  82,  p,  ffn2,  SNRlt  SNR2  ]. 

Six  derivatives  have  to  be  computed-  They  are 
3R/38J  »  OA/ae^SA8  +  AS(3Aa/381), 

3R/382  -  (3A/392)SAH  ♦  AS(3AH/382), 

3R/3p  .  A(3S/3p)AH, 

3R/3(ffn2)  -  A(3S/3(cn2)AB  ♦  IB, 

3R/3SNRJ  -  A( 3S/ 3SNRX )A8 
3R/3SNR2  -  A(3S/3SNR2)AH 

Note  that 


?S/2p  -  a2  io(SN8l-fSNR2)/20) 
•  « 

3S/3(«n2)  - 


0  1 
1  0 


3S/3SNRX  - 


as/asNR^  ■ 


■  10(SNR1)/10)  p  10(SNR1+SNR2)/20) 

p  10(SNR1*SNR2)/20)  10(SNR2)/10)  t 

((Ln(10)/10)  ioCSNR^/IO)  e((Ln(lC)/20)10<SNKKSNTR2>  20> 
[  p((Ln(10)/20)10^SNRl*SNK2>/20>  0 

0  D((Ln(10)  :o)io(SNRi‘SNR2)-:o> 

.  p((I^(10)/20,10<snr1*snr2>/20'  ((Ln<10)/10)  ic-< SKR1  >  10> 


1'9 


and 


OA/d^) 


0  0 

(coA/c)cos(©i)ej  ♦!  0 


[  (m-l)(a>A/c)cos( 8^)e^m-^#l  0  J 

r°  o  .  I 

0  (a>A/c)cos(  ©2)^  ^2 


(dA/302)  ~ 


L  0  (m-l)(wA/c)cos(e2)ej(®-1U2  J  . 

Having  defined  all  the  above  quantities,  it  is  easy  to  do  the  multiplica¬ 
tions  needed  in  equation  (A. 7)  and  determine  the  CRLB. 
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